{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 作业内容\n",
    "\n",
    "> 如果是0.1mm厚的金箔，6MeV的阿尔法粒子穿过时平均散射次数为1.7次。根据泊松分布（设期望为$\\lambda = 1.7$）抽取散射次数，并对每次散射，使用微分散射截面中的角度依赖$p(\\theta)=\\frac{\\sin \\theta}{\\sin^4{\\theta/2}}$抽样散射角，最后统计100万个阿尔法粒子穿过此金箔后沿束流方向的动量分布。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 参考答案"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import quad\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'times of scattering')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAHiCAYAAADF4pQuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAcRklEQVR4nO3de7CkdX3n8fdHBm+ggmEkOIBDDEnEVAXZCWLMptxguCbB1MYtiDHEmJ1sClOaa6GVWkwiu1hlNGXFWIVhFC8RiZdIhMSg0bhuVmVARAGNIw7OwAiTcI+JEfzuH/2bpMEzc86cOXN6vmfer6qu0/17nu7+PT2X9zxPP92TqkKSJO37HjXrCUiSpIUx2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbWiJJjk7yQJIDZj2XnUnyM0m2jHk+a9bzWawkf5Xk3FnPQ1pu8XPa0uIk2Qz8clV9ZNZzWagkXwF+o6o+uEzP9zZga1X97tTYZpq9btK+wj1taf/yNODGWU9isTLh31vab/mbX1qEJO8Ajgb+chxq/p0ka5NUklVjnY8neU2Svx/r/GWS70ryriT3Jbkmydqpx/yBJFcnuSvJl5L8t6llZyS5Kcn9SW5L8ls7mdejkvxukluT3Jnk7UmelOQxSR4ADgA+N/a4H3nfJHnDuN+9SW5I8oNj2eOS/OF43HuTfDLJ48ayP0/y9TH+iSTPHOPrgRcBvzO1/d/xuo11Txqv0z1JPpfkeVPz+niSC5P8X+AbwPeMsV8ey39xzOd1Se5O8tUkp0/d/5gxr/uTfCTJm5K8c7d/0aV9QVV58eJlERdgM/D8qdtrgQJWjdsfBzYBTweeBNwE/APwfGAV8HbgrWPdg4AtwEvGshOAfwSeOZZvA/7zuH4ocMJO5vRL4zm/BzgYeD/wjqnlBXzvTu57KnAtcAgQ4BnAEWPZm8b2rGES/h8BHjP1nE8AHgP8EXD91GO+DXjNPK/bGuCfgDOY7Ej8xLi9eup1/BrwzPHaHDjGfnks/0XgW8B/H3P7VeB2/uPtv/8HvA54NPCjwH3AO2f9+8eLl8Vc3NOW9q63VtVXqupe4K+Ar1TVR6rqQeDPgR0ng/0ksLmq3lpVD1bVdcD7gJ8dy78FHJfkiVV191g+lxcBr6+qW6rqAeCVwNk79v7n8S0m8f0BJsG7uaq2jcPRvwS8vKpuq6qHqurvq+qbAFW1oaruH7dfDfxQkiftxmv088BVVXVVVX27qq4GNjKJ+A5vq6obx2vzrTke49aqektVPQRcChwBHJ7kaOCHgf9ZVf9WVZ8ErtiNuUn7FKMt7V13TF3/lzluHzyuPw149jg8fE+Se5gE+LvH8v/KJGK3Jvm7JM/ZyfM9Fbh16vatTPZOD59volX1t8AfM9mrviPJxUmeCBwGPBaY65D6AUkuSvKVJPcx2Ytm3Gehnga88BHb/qNMwrvDlnke4+tT2/GNcfVgJq/HXVNjC3ksaZ9ltKXFW8qPXmwB/q6qDpm6HFxVvwpQVddU1VnAU4C/AC7fyePcziSCOxwNPMjD/7GwU1X1xqr6T0wORX8f8NtMDtP/K5PD/I/0c8BZTA75P4nJWwQwObwOc79GjxzbwuQQ/vS2H1RVF+3iPgu1DXhyksdPjR21yMeSZs5oS4t3B5P3jpfCh4DvS/LiJAeOyw8neUaSRyd5UZInjUPD9wEP7eRx3g38+jj56mDgfwHvGYfjd2k837OTHAj8M5NQP1RV3wY2AK9P8tSxd/2cJI9hcjj9m0zeg378eL5pc71Gjxx7J/BTSU4dj/3YJM9LcuR8c55PVd3K5FD7q8fr+Bzgp/b0caVZMdrS4v1v4HfHId05z+ZeqKq6HzgFOJvJ3vLXgdcyObkL4MXA5nEI+n8weR94LhuAdwCfAL7KJLy/tsBpPBF4C3A3k8Pq/8TkBC6A3wI+D1wD3DXm9igmJ9PdCtzG5ES7Tz3iMS9h8l78PUn+Yow97HWrqi1M9tZfBWxnsuf92yzd308vAp4ztuc1wHuY/ENDascvV5G0X0nyHuCLVXXBrOci7S73tCWtaOOw/9PHZ9hPY7JX/xfz3U/aFy3kYyCS1Nl3M/m8+ncBW4FfrarPznZK0uJ4eFySpCY8PC5JUhNGW5KkJvbp97QPO+ywWrt27aynIUnSsrn22mv/sapWz7Vsn4722rVr2bhx46ynIUnSskly686WeXhckqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1JkppYNesJ6OHWnn/lrKcwr80XnTnrKUjSfsk9bUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKamDfaSR6b5DNJPpfkxiS/N8aPSfLpJF9O8p4kjx7jjxm3N43la6ce65Vj/EtJTt1bGyVJ0kq0kD3tbwI/XlU/BBwPnJbkJOC1wBuq6ljgbuClY/2XAndX1fcCbxjrkeQ44GzgmcBpwJ8kOWApN0aSpJVs3mjXxAPj5oHjUsCPA+8d45cCLxjXzxq3GctPTpIxfllVfbOqvgpsAk5ckq2QJGk/sKD3tJMckOR64E7gauArwD1V9eBYZSuwZlxfA2wBGMvvBb5renyO+0iSpHksKNpV9VBVHQ8cyWTv+BlzrTZ+ZifLdjb+MEnWJ9mYZOP27dsXMj1JkvYLu3X2eFXdA3wcOAk4JMmO/9rzSOD2cX0rcBTAWP4k4K7p8TnuM/0cF1fVuqpat3r16t2ZniRJK9pCzh5fneSQcf1xwPOBm4GPAT87VjsX+OC4fsW4zVj+t1VVY/zscXb5McCxwGeWakMkSVrpVs2/CkcAl44zvR8FXF5VH0pyE3BZktcAnwUuGetfArwjySYme9hnA1TVjUkuB24CHgTOq6qHlnZzJElaueaNdlXdADxrjvFbmOPs76r6V+CFO3msC4ELd3+akiTJb0STJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpiVWznoD6WXv+lbOewi5tvujMWU9BkvYK97QlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNTFvtJMcleRjSW5OcmOSl4/xVye5Lcn143LG1H1emWRTki8lOXVq/LQxtinJ+XtnkyRJWplWLWCdB4HfrKrrkjwBuDbJ1WPZG6rqddMrJzkOOBt4JvBU4CNJvm8sfhPwE8BW4JokV1TVTUuxIZIkrXTzRruqtgHbxvX7k9wMrNnFXc4CLquqbwJfTbIJOHEs21RVtwAkuWysa7QlSVqA3XpPO8la4FnAp8fQy5LckGRDkkPH2Bpgy9Tdto6xnY0/8jnWJ9mYZOP27dt3Z3qSJK1oC452koOB9wGvqKr7gDcDTweOZ7In/oc7Vp3j7rWL8YcPVF1cVeuqat3q1asXOj1Jkla8hbynTZIDmQT7XVX1foCqumNq+VuAD42bW4Gjpu5+JHD7uL6zcUmSNI+FnD0e4BLg5qp6/dT4EVOr/QzwhXH9CuDsJI9JcgxwLPAZ4Brg2CTHJHk0k5PVrliazZAkaeVbyJ72c4EXA59Pcv0YexVwTpLjmRzi3gz8CkBV3ZjkciYnmD0InFdVDwEkeRnwYeAAYENV3biE2yJJ0oq2kLPHP8nc70dftYv7XAhcOMf4Vbu6nyRJ2jm/EU2SpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJuaNdpKjknwsyc1Jbkzy8jH+5CRXJ/ny+HnoGE+SNybZlOSGJCdMPda5Y/0vJzl3722WJEkrz0L2tB8EfrOqngGcBJyX5DjgfOCjVXUs8NFxG+B04NhxWQ+8GSaRBy4Ang2cCFywI/SSJGl+80a7qrZV1XXj+v3AzcAa4Czg0rHapcALxvWzgLfXxKeAQ5IcAZwKXF1Vd1XV3cDVwGlLujWSJK1gu/WedpK1wLOATwOHV9U2mIQdeMpYbQ2wZepuW8fYzsYlSdICLDjaSQ4G3ge8oqru29Wqc4zVLsYf+Tzrk2xMsnH79u0LnZ4kSSvegqKd5EAmwX5XVb1/DN8xDnszft45xrcCR03d/Ujg9l2MP0xVXVxV66pq3erVq3dnWyRJWtEWcvZ4gEuAm6vq9VOLrgB2nAF+LvDBqfFfGGeRnwTcOw6ffxg4Jcmh4wS0U8aYJElagFULWOe5wIuBzye5foy9CrgIuDzJS4GvAS8cy64CzgA2Ad8AXgJQVXcl+QPgmrHe71fVXUuyFZIk7QfmjXZVfZK5348GOHmO9Qs4byePtQHYsDsTlCRJE34jmiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1MSqWU9gua09/8pZT0GSpEVxT1uSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNTFvtJNsSHJnki9Mjb06yW1Jrh+XM6aWvTLJpiRfSnLq1PhpY2xTkvOXflMkSVrZFrKn/TbgtDnG31BVx4/LVQBJjgPOBp457vMnSQ5IcgDwJuB04DjgnLGuJElaoFXzrVBVn0iydoGPdxZwWVV9E/hqkk3AiWPZpqq6BSDJZWPdm3Z7xpIk7af25D3tlyW5YRw+P3SMrQG2TK2zdYztbFySJC3QYqP9ZuDpwPHANuAPx3jmWLd2Mf4dkqxPsjHJxu3bty9yepIkrTyLinZV3VFVD1XVt4G38B+HwLcCR02teiRw+y7G53rsi6tqXVWtW7169WKmJ0nSijTve9pzSXJEVW0bN38G2HFm+RXAnyV5PfBU4FjgM0z2tI9NcgxwG5OT1X5uTyYu7cza86+c9RTmtfmiM2c9BUkNzRvtJO8GngcclmQrcAHwvCTHMznEvRn4FYCqujHJ5UxOMHsQOK+qHhqP8zLgw8ABwIaqunHJt0aSpBVsIWePnzPH8CW7WP9C4MI5xq8Crtqt2UmSpH/nN6JJktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpiXmjnWRDkjuTfGFq7MlJrk7y5fHz0DGeJG9MsinJDUlOmLrPuWP9Lyc5d+9sjiRJK9dC9rTfBpz2iLHzgY9W1bHAR8dtgNOBY8dlPfBmmEQeuAB4NnAicMGO0EuSpIWZN9pV9QngrkcMnwVcOq5fCrxgavztNfEp4JAkRwCnAldX1V1VdTdwNd/5DwFJkrQLi31P+/Cq2gYwfj5ljK8Btkytt3WM7Wz8OyRZn2Rjko3bt29f5PQkSVp5lvpEtMwxVrsY/87Bqoural1VrVu9evWSTk6SpM4WG+07xmFvxs87x/hW4Kip9Y4Ebt/FuCRJWqDFRvsKYMcZ4OcCH5wa/4VxFvlJwL3j8PmHgVOSHDpOQDtljEmSpAVaNd8KSd4NPA84LMlWJmeBXwRcnuSlwNeAF47VrwLOADYB3wBeAlBVdyX5A+Casd7vV9UjT26TJEm7MG+0q+qcnSw6eY51CzhvJ4+zAdiwW7OTJEn/zm9EkySpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU2smvUEJO2b1p5/5aynsEubLzpz1lOQlp172pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpiT2KdpLNST6f5PokG8fYk5NcneTL4+ehYzxJ3phkU5IbkpywFBsgSdL+Yin2tP9LVR1fVevG7fOBj1bVscBHx22A04Fjx2U98OYleG5JkvYbe+Pw+FnApeP6pcALpsbfXhOfAg5JcsReeH5JklakPY12AX+T5Nok68fY4VW1DWD8fMoYXwNsmbrv1jEmSZIWYNUe3v+5VXV7kqcAVyf54i7WzRxj9R0rTeK/HuDoo4/ew+lJkrRy7NGedlXdPn7eCXwAOBG4Y8dh7/HzzrH6VuCoqbsfCdw+x2NeXFXrqmrd6tWr92R6kiStKIuOdpKDkjxhx3XgFOALwBXAuWO1c4EPjutXAL8wziI/Cbh3x2F0SZI0vz05PH448IEkOx7nz6rqr5NcA1ye5KXA14AXjvWvAs4ANgHfAF6yB88tSdJ+Z9HRrqpbgB+aY/yfgJPnGC/gvMU+nyRJ+zu/EU2SpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmjLYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKmJVbOegCQtxtrzr5z1FOa1+aIzZz0FrTDuaUuS1ITRliSpCaMtSVITRluSpCaMtiRJTRhtSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJoy2JElNGG1Jkpow2pIkNWG0JUlqwmhLktSE0ZYkqQmjLUlSE0ZbkqQmVs16ApK0Uq09/8pZT2GXNl905qynoN3knrYkSU0YbUmSmjDakiQ1YbQlSWrCaEuS1MSyRzvJaUm+lGRTkvOX+/klSepqWaOd5ADgTcDpwHHAOUmOW845SJLU1XJ/TvtEYFNV3QKQ5DLgLOCmZZ6HJO339vXPkYOfJX+k5T48vgbYMnV76xiTJEnzWO497cwxVg9bIVkPrB83H0jypSWew2HAPy7xY87CStkO2A+3Ja9dhpnsuZXy67JStgP2w23ZT/+sPG1nC5Y72luBo6ZuHwncPr1CVV0MXLy3JpBkY1Wt21uPv1xWynaA27KvWinbslK2A9yWfdVybstyHx6/Bjg2yTFJHg2cDVyxzHOQJKmlZd3TrqoHk7wM+DBwALChqm5czjlIktTVsv8vX1V1FXDVcj/vlL126H2ZrZTtALdlX7VStmWlbAe4LfuqZduWVNX8a0mSpJnza0wlSWpiv4n2Svn61CQbktyZ5AuznsueSnJUko8luTnJjUlePus5LVaSxyb5TJLPjW35vVnPaU8kOSDJZ5N8aNZz2RNJNif5fJLrk2yc9Xz2RJJDkrw3yRfHn5nnzHpOuyvJ949fix2X+5K8YtbzWqwkvz7+vH8hybuTPHavP+f+cHh8fH3qPwA/weRjZ9cA51RVu29iS/JjwAPA26vqB2c9nz2R5AjgiKq6LskTgGuBFzT9dQlwUFU9kORA4JPAy6vqUzOe2qIk+Q1gHfDEqvrJWc9nsZJsBtZVVfvPNie5FPg/VfWn49M3j6+qe2Y9r8Uafy/fBjy7qm6d9Xx2V5I1TP6cH1dV/5LkcuCqqnrb3nze/WVP+9+/PrWq/g3Y8fWp7VTVJ4C7Zj2PpVBV26rqunH9fuBmmn5DXk08MG4eOC4t/0Wc5EjgTOBPZz0XTSR5IvBjwCUAVfVvnYM9nAx8pWOwp6wCHpdkFfB4HvG9I3vD/hJtvz51H5dkLfAs4NOzncnijUPK1wN3AldXVddt+SPgd4Bvz3oiS6CAv0ly7fi2xa6+B9gOvHW8bfGnSQ6a9aT20NnAu2c9icWqqtuA1wFfA7YB91bV3+zt591foj3v16dqdpIcDLwPeEVV3Tfr+SxWVT1UVccz+aa/E5O0e/siyU8Cd1bVtbOeyxJ5blWdwOR/FjxvvL3U0SrgBODNVfUs4J+BzufmPBr4aeDPZz2XxUpyKJMjtscATwUOSvLze/t595doz/v1qZqN8f7v+4B3VdX7Zz2fpTAOW34cOG3GU1mM5wI/Pd4Lvgz48STvnO2UFq+qbh8/7wQ+wOStso62Alunjt68l0nEuzoduK6q7pj1RPbA84GvVtX2qvoW8H7gR/b2k+4v0fbrU/dB4+StS4Cbq+r1s57PnkiyOskh4/rjmPyB/uJsZ7X7quqVVXVkVa1l8ufkb6tqr+897A1JDhonODIOJZ8CtPzURVV9HdiS5PvH0Mn0/i+Nz6HxofHha8BJSR4//i47mcl5OXvVsn8j2iyspK9PTfJu4HnAYUm2AhdU1SWzndWiPRd4MfD58V4wwKvGt+Z1cwRw6Tgj9lHA5VXV+uNSK8DhwAcmf5+yCvizqvrr2U5pj/wa8K6x43EL8JIZz2dRkjyeySd5fmXWc9kTVfXpJO8FrgMeBD7LMnwz2n7xkS9JklaC/eXwuCRJ7RltSZKaMNqSJDVhtCVJasJoS5LUhNGWJKkJoy1JUhNGW5KkJv4/qkmDEOmEDf0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#sampling times of scattering\n",
    "n = 10000\n",
    "times = np.random.poisson(1.7, (n,))\n",
    "\n",
    "#check\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(8, 8)\n",
    "ax.hist(times)\n",
    "ax.set_title('times of scattering')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "化简$p(\\theta)=\\frac{\\sin \\theta}{\\sin^4{\\theta/2}}$，显然有$p(\\theta)=\\frac{2\\sin(\\theta/2)\\cos(\\theta/2)}{\\sin^4(\\theta/2)}$，易得：\n",
    "\n",
    "$$p(\\theta)d\\theta=4\\frac{d\\sin (\\theta/2)}{\\sin^3(\\theta/2)}$$\n",
    "\n",
    "令$x = \\sin(\\theta/2)$，那么就可以对$p(x)=Ax^{-3}$抽样"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 蒙特卡洛舍选法\n",
    "\n",
    "遇到不会做的可以无脑选这个，反正总是能抽出来的"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#monte-carlo sampling function\n",
    "def sampling(func, range_=[[0.1, 1], [0, 1]], n=1000):\n",
    "    r'''\n",
    "    Monte-Carlo sampling.\n",
    "    Args: 3\n",
    "        func: a function, the pdf. to be sampled.\n",
    "        range_: a list, [[xmin, xmax], [ymin, ymax]].\n",
    "        n: an integer, how many variables you want to simulate.\n",
    "    Return: 1\n",
    "        res: an array with n variables respect to the pdf. func.\n",
    "    '''\n",
    "    [xmin, xmax], [ymin, ymax] = range_\n",
    "    samples = np.zeros(n, dtype=np.float32)\n",
    "    i = 0\n",
    "    while i < n:\n",
    "        x = np.random.uniform(xmin, xmax)\n",
    "        y = np.random.uniform(ymin, ymax)\n",
    "        if y < func(x):\n",
    "            samples[i] = x\n",
    "            i = i + 1\n",
    "    return samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Amplitude: 49.499999999999986\n",
      "Maximum of p(k): 20.202020202020204\n",
      "Minimum of p(k): 0.020202020202020207\n"
     ]
    }
   ],
   "source": [
    "#define the pdf. function\n",
    "def raw_func(x):\n",
    "    return x**(-3)\n",
    "\n",
    "amplitude, _ = quad(raw_func, 0.1, 1)\n",
    "\n",
    "def func(x):\n",
    "    global amplitude\n",
    "    return raw_func(x) / amplitude\n",
    "\n",
    "#check\n",
    "x0 = np.linspace(0.1, 1, 100)\n",
    "px0 = func(x0)\n",
    "pmax = px0.max()\n",
    "pmin = px0.min()\n",
    "print('Amplitude: {}'.format(amplitude))\n",
    "print('Maximum of p(k): {}'.format(pmax))\n",
    "print('Minimum of p(k): {}'.format(pmin))\n",
    "\n",
    "range_ = [[0.1, 1], [pmin, pmax]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAHjCAYAAADsYlzyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZze873//8drlkwy2XeSIEGEIIJYo8RSJfFFtbailpI6h9OWtqfpckpPT91UnbYcPfWjlhbFadVWoUKpIqpJRUSCoCERJBKyTSbbvH9/XNeMScwkk9mu65rrcb/dPrfrc33W13yMPOf9Wd6fSCkhSZLyU0muC5AkSY0zqCVJymMGtSRJecygliQpjxnUkiTlMYNakqQ8ZlBLkpTHDGpJkvKYQS11cBHxckSMy3UdkprHoJY6uJTS7imlJ5u6fETsHRHPRERVRDwfEdvXmzcoIhZEREVE3BQRb0XEioh4ISKObZMfQCpyBrWkOhExBJgM/BjoC7wJfK/eIuOBR4AyYD5wGNAT+A/g/yJiaDuWKxUFg1rqICLiWxHxTraF+2pEHJmdPi8ijqo3/o2ImBkRyyLi7ojoXG8z/w3cmFJ6IKW0GrgL2K/e/PHA5JTSqpTS5SmleSmlmpTSH4F/Avu2z08rFQ+DWuoAImIEcDGwX0qpO/AZYF4ji58CHAMMA0YB52S30QM4AfhVvWVLgOrs/HLgUGBKA/sfCOwCvNziH0bSRspyXYCkVrEBqABGRsTilNK8zSx7bUppIUBEPAiMzk4/EigHZkZE7bIVwP3Z8UOBF1NKK+pvLBvgdwC/Tim90go/i6R6bFFLHUBK6XXga8DlwKKIuCsiBjWy+Hv1xquAbtnxocADKaVetQPwBJlr0pA97V1/QxFRAtwGrCXTopfUygxqqYNIKf02pXQIsAOQyNwQtjUqyAQ3ABExDBgDPJCdNB54qN78AG4CBgKfSymta371khpjUEsdQESMiIgjIqKCzDXl1WROh2+NvwOHZR/B2g74LfDdlNLSbGhXbHJq+5fAbsD/y954JqkNGNRSx1ABXAl8QObU9gDgO1u5jT8DDwKvAU8Dt6WUbszOm0C9094RsQPwZTLXt9+LiJXZ4YwW/RSSPiFSSrmuQVKei4jJwHUppclbXFhSq7JFLakpniRzY5mkdmaLWpKkPGaLWpKkPGZQS5KUx/KyZ7J+/fqloUOH5roMSZLaxfTp0z9IKfVvaF5eBvXQoUOZNm1arsuQJKldRMRbjc3z1LckSXnMoJYkKY8Z1JIk5bG8vEYtSbXWrVvHggULqK6uznUpUot17tyZIUOGUF5e3uR1DGpJeW3BggV0796doUOHUu892VLBSSmxZMkSFixYwLBhw5q8nqe+JeW16upq+vbta0ir4EUEffv23eqzQwa1pLxnSKujaM7vskEtSVIe8xq1JG3GkiVLOPLIIwF47733KC0tpX///sybN49BgwYxe/bsHFeojs4WtSRtRt++fZkxYwYzZszgwgsv5JJLLqn7XlLS+v+Erl+/vtW3qcJmi1pS4fja12DGjNbd5ujR8POfN2vVDRs2cMEFF/Dss88yePBg7r//frp06cIbb7zBRRddxOLFi6msrOTGG29k11135a233uK8885j8eLF9O/fn1tuuYXtt9+ec845hz59+vDCCy8wevRo/vjHP/Lss8/Sv39/ampq2GWXXXjuuefo169f6/7sKgi2qCWpmebOnctFF13Eyy+/TK9evbjnnnsAmDhxIv/zP//D9OnTufrqq/nXf/1XAC6++GK++MUvMnPmTM444wy+8pWv1G3rtdde47HHHuNnP/sZZ555JnfccQcAjz32GHvttZchXcRsUUsqHM1s+baVYcOGMXr0aAD23Xdf5s2bx8qVK3n22Wc5+eST65Zbs2YNAFOnTuUPf/gDAGeddRb//u//XrfMySefTGlpKQDnnXceJ5xwAl/72te4+eabOffcc9vrR1IeMqglqZkqKirqxktLS1m9ejU1NTX06tWLGU04RV//UZ2uXbvWjW+33XYMHDiQP//5z/ztb3+ra12rOHnqW5JaUY8ePRg2bBi/+93vgExvVC+++CIABx98MHfddRcAd9xxB4ccckij2zn//PM588wzOeWUU+pa2ipOBrUktbI77riDm266ib322ovdd9+d+++/H4Brr72WW265hVGjRnHbbbdxzTXXNLqN448/npUrV3raW0RKKdc1fMKYMWPStGnTcl2GpDwwZ84cdtttt1yX0e6mTZvGJZdcwl//+tdcl6JW1tDvdERMTymNaWj5jn+NeuFCKCuDAQNyXYkkNcmVV17JL3/5S69NCyiGU9/77w/f+U6uq5CkJps0aRJvvfXWZq9hq3h0/KDu2hVWrcp1FZIkNcsWgzoitouIJyJiTkS8HBFfzU7vExFTImJu9rN3I+ufnV1mbkSc3do/wBYZ1JKkAtaUFvV64Osppd2AA4GLImIkMAl4PKU0HHg8+30jEdEHuAw4ANgfuKyxQG8zBrUkqYBtMahTSu+mlP6RHV8BzAEGAycAv84u9mvgxAZW/wwwJaW0NKX0ITAFOKY1Cm8yg1qSVMC26q7viBgK7A38DRiYUnoXMmEeEQ3dVj0YmF/v+4LstIa2PRGYCLD99ttvTVmb17UrvPNO621PUk4NnfRQq25v3pUTtrhMRHDmmWdy2223AZk3XG277bYccMAB/PGPf9zqfX700Uf89re/resDfGs8/PDD/Md//AerVq0ipcRxxx3H1Vdf3eT1n3zySa6++uom111VVcUFF1zAzJkzSSnRq1cvHnnkEbp167bVtTdVt27dWLlyJQsXLuQrX/kKv//971u8zVtvvZVvfvObDB48mOrqar785S9zySWXtEK1Hxs3bhxXX301Y8aM4YorruA7rXQjc5NvJouIbsA9wNdSSsubuloD0xp8cDuldENKaUxKaUz//v2bWtaW2aKW1EJdu3Zl1qxZrF69GoApU6YweHCDbY4m+eijj/jf//3frV5v1qxZXHzxxdx+++3MmTOHWbNmseOOOzZ5/ea8QvOaa65h4MCBvPTSS8yaNYubbrqJ8vLyrd5OcwwaNKhVQrrWqaeeyowZM3jmmWf40Y9+xPz587e8UjNdccUVrbatJgV1RJSTCek7Ukp/yE5+PyK2zc7fFljUwKoLgO3qfR8CLGx+uc1gUEtqBcceeywPPZRpzd95552cfvrpdfOWLl3KiSeeyKhRozjwwAOZOXMmAJdffjnnnXce48aNY8cdd+Taa68FMo9fvfHGG4wePZpvfvObAPzkJz9hv/32Y9SoUVx22WUN1nDVVVfx3e9+l1133RWAsrKyulb5gw8+yAEHHMDee+/NUUcdxfvvv19Xw8SJEzn66KP54he/uNH2Gqu7vnfffXejP0pGjBhR18f5iSeeyL777svuu+/ODTfcULdMt27d+Na3vsW+++7LUUcdxfPPP193DB544AEg08I94YQTOOaYYxgxYgQ/+MEPPrHvefPmsccee9Qtf9JJJ3HMMccwfPjwjV5octNNN7HLLrswbtw4LrjgAi6++OIGj1+tvn37svPOO/Puu+8CsHjxYj73uc+x3377sd9++/HMM88A8Je//IXRo0czevRo9t57b1asWMGTTz7JcccdV7etiy++mFtvvXWj7U+aNInVq1czevRozjjjjM3W0hRNues7gJuAOSmln9ab9QBQexf32cD9Daz+J+DoiOidvYns6Oy09mNQS2oFp512GnfddRfV1dXMnDmTAw44oG7eZZddxt57783MmTO54oorNgrEV155hT/96U88//zz/OAHP2DdunVceeWV7LTTTsyYMYOf/OQnPProo8ydO5fnn3+eGTNmMH36dJ566qlP1DBr1iz23XffBus75JBDeO6553jhhRc47bTTuOqqq+rmTZ8+nfvvv5/f/va3G62zubprnXfeefz4xz/moIMO4nvf+x5z586tm3fzzTczffp0pk2bxrXXXsuSJUsAWLVqFePGjWP69Ol0796d733ve0yZMoV7772X73//+3XrP//889xxxx3MmDGD3/3ud2ypR8oZM2Zw991389JLL3H33Xczf/58Fi5cyA9/+EOee+45pkyZwiuvvLLZbQC8/fbbVFdXM2rUKAC++tWvcskll/D3v/+de+65h/PPPx+Aq6++ml/84hfMmDGDv/71r3Tp0mWL24ZMhzVdunRhxowZrdJpTVOuUY8FzgJeioja18F8B7gS+L+I+BLwNnAyQESMAS5MKZ2fUloaET8E/p5d7z9TSktbXPXW6NoVqqqgpgZKOv5j45LaxqhRo5g3bx533nkn48eP32je008/Xfcu6iOOOIIlS5awbNkyACZMmEBFRQUVFRUMGDCgrqVb36OPPsqjjz7K3nvvDcDKlSuZO3cuhx56aJPrW7BgAaeeeirvvvsua9euZdiwYXXzjj/++AZDprG6e/bsWbfM6NGjefPNN3n00Ud57LHH2G+//Zg6dSq77bYb1157Lffeey8A8+fPZ+7cufTt25dOnTpxzDGZ+4b33HNPKioqKC8vZ88992TevHl12/70pz9N3759ATjppJN4+umnGTOmwV40ATjyyCPrahs5ciRvvfUWH3zwAYcddhh9+vQBMq8Lfe211xpc/+677+aJJ57g1Vdf5cYbb6Rz585A5p3fs2fPrltu+fLlrFixgrFjx3LppZdyxhlncNJJJzFkyJBGa2tLWwzqlNLTNHytGeDIBpafBpxf7/vNwM3NLbDFunaFlGD16sy4JDXT8ccfzze+8Q2efPLJutYjZN6QtanaV1hu+irMhq4Tp5T49re/zZe//OWNpv/iF7/gxhtvBGDy5MnsvvvuTJ8+nb322usT2/i3f/s3Lr30Uo4//niefPJJLr/88rp5XRv5t29zddfXrVs3TjrpJE466SRKSkqYPHky77//Po899hhTp06lsrKScePGUV1dDUB5eXnddkpKSuqOQUlJyUY//6b7amjf9TV0LLfmfRWnnnoq1113HVOnTmXChAkce+yxbLPNNtTU1DB16tRP/DEzadIkJkyYwOTJkznwwAN57LHHKCsro6ampm6Z2p+5LXX8JmbtL6invyW10Hnnncf3v/999txzz42mH3rooXWnOJ988kn69etHjx49Gt1O9+7dWbFiRd33z3zmM9x8882sXLkSgHfeeYdFixZx0UUXMWPGDGbMmMGgQYP45je/yRVXXFHXYqypqeGnP81ckVy2bFndteRf//rXNEVT6n7mmWf48MMPAVi7di2zZ89mhx12YNmyZfTu3ZvKykpeeeUVnnvuuSbts74pU6awdOlSVq9ezX333cfYsWO3ehv7778/f/nLX/jwww9Zv3593RmCzTnooIM466yz6t5edvTRR3PdddfVza99l/gbb7zBnnvuybe+9S3GjBnDK6+8wg477MDs2bNZs2YNy5Yt4/HHH29wH+Xl5axbt26rf56GdPyXchjUUofSlMep2sqQIUP46le/+onpl19+Oeeeey6jRo2isrJyi0HZt29fxo4dyx577MGxxx7LT37yE+bMmcNBBx0EZFqwt99+OwM2eZnQqFGj+PnPf87pp59OVVUVEcGECRPqajj55JMZPHgwBx54IP/85z+3+PM0pe433niDf/mXfyGlRE1NDRMmTOBzn/sca9eu5frrr2fUqFGMGDGCAw88cIv729QhhxzCWWedxeuvv84XvvCFzZ72bszgwYP5zne+wwEHHMCgQYMYOXLkRqfuG/Otb32LffbZh+985ztce+21XHTRRYwaNYr169dz6KGHcv311/Pzn/+cJ554gtLSUkaOHMmxxx5LRUUFp5xyCqNGjWL48OF1lys2NXHiREaNGsU+++zT4uvUHf81l7/7HZxyCrz0EmTvHpRUOIr1NZcd3a233sq0adM2ask218qVK+nWrRvr16/ns5/9LOeddx6f/exnW6HKtrG1r7n01LckqaBdfvnljB49mj322INhw4Zx4okNdZRZuDz1LUlqd+eccw7nnHNOq2xra3pmK0S2qCXlvXy8RCc1R3N+lw1qSXmtc+fOLFmyxLBWwUspsWTJkrrnt5vKU9+S8tqQIUNYsGABixcvznUpUot17tx5qztOMagl5bXy8vKNetmSio2nviVJymMdP6grKjJ9fBvUkqQC1PGDOsI3aEmSClbHD2owqCVJBcugliQpjxnUkiTlMYNakqQ8ZlBLkpTHDGpJkvKYQS1JUh4zqCVJymMGtSRJeay4gtrX5EmSCkzxBHVNDaxZk+tKJEnaKsUT1ODpb0lSwTGoJUnKYwa1JEl5zKCWJCmPGdSSJOUxg1qSpDxmUEuSlMeKK6hXrsxtHZIkbaXiCmpb1JKkAlMUQb3blX8F4MrfT2fopIdyXI0kSU1XFEFdXd4JgC7rqnNciSRJW6cogjpFCVXlFVQa1JKkAlMUQQ1QVd7ZoJYkFZyiCerV5Z3pss63Z0mSCkvRBPUqW9SSpAJUNEG9urwzlWsNaklSYSmaoK7qVOGpb0lSwSmeoPbUtySpABVNUGduJjOoJUmFpWxLC0TEzcBxwKKU0h7ZaXcDI7KL9AI+SimNbmDdecAKYAOwPqU0ppXq3mq2qCVJhWiLQQ3cClwH/KZ2Qkrp1NrxiPhvYNlm1j88pfRBcwtsLavLK6j0GrUkqcBsMahTSk9FxNCG5kVEAKcAR7RuWa2vylPfkqQC1NJr1J8C3k8pzW1kfgIejYjpETFxcxuKiIkRMS0ipi1evLiFZX1SVXlnyms2UL5hXatvW5KkttLSoD4duHMz88emlPYBjgUuiohDG1swpXRDSmlMSmlM//79W1jWJ60u7wzgI1qSpILS7KCOiDLgJODuxpZJKS3Mfi4C7gX2b+7+WqqqvALATk8kSQWlJS3qo4BXUkoLGpoZEV0jonvtOHA0MKsF+2uRqk6ZFrV3fkuSCskWgzoi7gSmAiMiYkFEfCk76zQ2Oe0dEYMiYnL260Dg6Yh4EXgeeCil9Ejrlb51Pj71bVBLkgpHU+76Pr2R6ec0MG0hMD47/iawVwvrazVV5baoJUmFp4h6Jsteo/ZmMklSASmaoK7y1LckqQAVXVB76luSVEiKJqhX1wW1p74lSYWjaIK69jnqLj5HLUkqIEUT1B/fTGZQS5IKR9EEdU1JKdVlnbyZTJJUUIomqKH2ndReo5YkFY4iC+oKT31LkgpKUQX16vLOVK5dnesyJElqsqIKak99S5IKTVEF9eryCm8mkyQVlKIK6kyL2qCWJBWOIgxqT31LkgpHUQX16vLOnvqWJBWUogrqqk4+niVJKixFFdS2qCVJhaaogrqqvDMVG9bD+vW5LkWSpCYpsqDOvJiDVatyW4gkSU1UVEFd+05qg1qSVCiKKqirDGpJUoExqCVJymNluS6gPa3OXqP+3NVTmD7knbrp866ckKuSJEnarOJqUXfKtKh9llqSVCiKKqhrbyYzqCVJhaKogrr2GnUX+/uWJBWIIgvqzDVqW9SSpEJRVEFde+q7y1qDWpJUGIoqqKu8Ri1JKjBFFdTrS8tYW1LmO6klSQWjqIIaMs9S+wYtSVKhKLqgrirv7KlvSVLBKL6g7tTZU9+SpIJRfEFd3pnKtatzXYYkSU1SdEG9sqKSrga1JKlAFF1Qr6joSo81vj1LklQYii6ol1d0pbtBLUkqEEUX1CsqKulRbVBLkgpD0QX18oqudFu7mkg1uS5FkqQtKsqgLiHRfU1VrkuRJGmLii6oV1R0BTCoJUkFoeiCennn2qD2OrUkKf9tMagj4uaIWBQRs+pNuzwi3omIGdlhfCPrHhMRr0bE6xExqTULb67aFrWPaEmSCkFTWtS3Asc0MP1nKaXR2WHypjMjohT4BXAsMBI4PSJGtqTY1rCiohKwRS1JKgxbDOqU0lPA0mZse3/g9ZTSmymltcBdwAnN2E6rWl7bovYRLUlSAWjJNeqLI2Jm9tR47wbmDwbm1/u+IDutQRExMSKmRcS0xYsXt6Cszfv4ZjKDWpKU/5ob1L8EdgJGA+8C/93AMtHAtNTYBlNKN6SUxqSUxvTv37+ZZW2Zd31LkgpJs4I6pfR+SmlDSqkGuJHMae5NLQC2q/d9CLCwOftrTWvLyqku6+TNZJKkgtCsoI6Ibet9/Swwq4HF/g4Mj4hhEdEJOA14oDn7a20rKio99S1JKghlW1ogIu4ExgH9ImIBcBkwLiJGkzmVPQ/4cnbZQcCvUkrjU0rrI+Ji4E9AKXBzSunlNvkpttLyim7eTCZJKghbDOqU0ukNTL6pkWUXAuPrfZ8MfOLRrVxbUVHpqW9JUkEoup7JoPZVl95MJknKf0Ub1LaoJUmFoDiDunNXbyaTJBWEogzqFbaoJUkFokiDupLO69dSvmFdrkuRJGmzijKol9s7mSSpQBRlUNe96rJ6ZY4rkSRp84oyqJd3tkUtSSoMRRnUdS1qbyiTJOW5Ig3qSsBXXUqS8l9RBvXyim6AQS1Jyn9FGdS1LWpfzCFJyndFG9Q1BD28mUySlOeKMqhTlLCyUxd6rPHxLElSfivKoIba/r5tUUuS8lvRBvWKCl/MIUnKf0Uc1JU+Ry1JyntFG9TLKzz1LUnKf0Ub1CsqutrXtyQp7xVtUHszmSSpEBRtUNfdTJZSrkuRJKlRRRzUlZSlGirXVee6FEmSGlW0QW1/35KkQlC0QW1/35KkQlC0Qb08+05qbyiTJOWzog3qFdmgtr9vSVI+K9qgXt65NqhtUUuS8lfxBnXdqW+vUUuS8lfRBnXdzWQGtSQpjxVtUFeXVbC2pMwWtSQprxVtUBOReYOWj2dJkvJY8QY19vctScp/RR3Udf19S5KUp4o8qCu9mUySlNeKOqiXV3SzRS1JymtFHdTeTCZJyndFHdTLK7rSfa03k0mS8ldRB/WKiq50W7sa1q/PdSmSJDWo6IMagOXLc1uIJEmNKOqgrn0xB8uW5bYQSZIaUdRBXdvfNx99lNtCJElqRFEH9fKKbpkRW9SSpDy1xaCOiJsjYlFEzKo37ScR8UpEzIyIeyOiVyPrzouIlyJiRkRMa83CW8Py2ha1QS1JylNNaVHfChyzybQpwB4ppVHAa8C3N7P+4Sml0SmlMc0rse0s75xtUXvqW5KUp7YY1Cmlp4Clm0x7NKVU+0zTc8CQNqitza2wRS1JynOtcY36PODhRuYl4NGImB4RE1thX61qZSdvJpMk5beylqwcEd8F1gN3NLLI2JTSwogYAEyJiFeyLfSGtjURmAiw/fbbt6SsJltfWsaq8s50tUUtScpTzW5RR8TZwHHAGSml1NAyKaWF2c9FwL3A/o1tL6V0Q0ppTEppTP/+/Ztb1lZbUVHpqW9JUt5qVlBHxDHAt4DjU0oNdpYdEV0jonvtOHA0MKuhZXNpeUU3T31LkvJWUx7PuhOYCoyIiAUR8SXgOqA7mdPZMyLi+uyygyJicnbVgcDTEfEi8DzwUErpkTb5KVrAFrUkKZ9t8Rp1Sun0Bibf1MiyC4Hx2fE3gb1aVF07WN65qy1qSVLeKuqeyQA+6twdlizJdRmSJDWo6IN6aWVPWLw412VIktSgog/qJZU9YeVKqK7OdSmSJH1C0Qf10i49MiO2qiVJecigruyZGTGoJUl5qOiDeolBLUnKY0Uf1LaoJUn5rOiD2ha1JCmfFX1QL6/oCmVlBrUkKS8VfVATAf36GdSSpLxkUAP0729QS5LykkENBrUkKW8Z1GBQS5LylkENBrUkKW8Z1JAJ6o8+gnXrcl2JJEkbMaghE9QAH3yQ2zokSdqEQQ0fB7WnvyVJecagBoNakpS3DGowqCVJecugBoNakpS3DGqAPn0yXYka1JKkPGNQA5SWQt++BrUkKe8Y1LV8MYckKQ8Z1LXsnUySlIcM6loGtSQpDxnUtQxqSVIeMqhr9e8PS5fChg25rkSSpDoGda3+/aGmBj78MNeVSJJUx6CuZacnkqQ8ZFDXMqglSXnIoK5lUEuS8pBBXcugliTlIYO6Vr9+mU+DWpKURwzqWp06Qc+eBrUkKa8Y1PXZ6YkkKc8Y1PUZ1JKkPGNQ12dQS5LyTFmuC8gHQyc9BMCV86o5/M0FDMxxPZIk1bJFXc/Syh70rloOKeW6FEmSAIN6I0u69KRTzXpYtizXpUiSBBjUG1la2TMz4nVqSVKeMKjrMaglSfnGoK5niUEtScozTQrqiLg5IhZFxKx60/pExJSImJv97N3Iumdnl5kbEWe3VuFtYWllj8yIQS1JyhNNbVHfChyzybRJwOMppeHA49nvG4mIPsBlwAHA/sBljQV6PljSxRa1JCm/NCmoU0pPAUs3mXwC8Ovs+K+BExtY9TPAlJTS0pTSh8AUPhn4eWNNeQWryjsb1JKkvNGSa9QDU0rvAmQ/BzSwzGBgfr3vC7LTPiEiJkbEtIiYtjiHQbm00hdzSJLyR1vfTBYNTGuwN5GU0g0ppTEppTH9a98NnQNLKnsY1JKkvNGSoH4/IrYFyH4uamCZBcB29b4PARa2YJ9tbmkXW9SSpPzRkqB+AKi9i/ts4P4GlvkTcHRE9M7eRHZ0dlre8tS3JCmfNPXxrDuBqcCIiFgQEV8CrgQ+HRFzgU9nvxMRYyLiVwAppaXAD4G/Z4f/zE7LW4u69Yb33oOamlyXIklS096elVI6vZFZRzaw7DTg/HrfbwZublZ1OfBu936wbl2mVT3Q92hJknLLnsk28V73fpmRBQtyW4gkSRjUn/CuQS1JyiMG9Sbe6943M2JQS5LygEG9iSWVPaG83KCWJOUFg3oTKUpg8GCDWpKUFwzqhgwZYlBLkvKCQd0Qg1qSlCcM6obUBnVqsFtySZLajUHdkCFDoLoaluZ1J2qSpCJgUDdkcPZNnJ7+liTlmEHdkCFDMp8GtSQpxwzqhtQG9Tvv5LYOSVLRM6gbss02UFJii1qSlHMGdUPKymDbbQ1qSVLOGdSN8VlqSVIeMKgbY1BLkvKAQd0Yg1qSlAcM6sYMGQIrVsDy5bmuRJJUxAzqxvgstSQpDxjUjTGoJUl5wKBujEEtScoDBnVjBg3KfBrUkqQcMqgb06fiH1kAABoGSURBVKkTDBxoUEuScsqg3hwf0ZIk5ZhBvTkGtSQpxwzqzRk82KCWJOWUQb05Q4bAhx/CqlW5rkSSVKQM6s3xvdSSpBwzqDfHZ6klSTlmUG+OLWpJUo4Z1JszeHDm0xa1JClHDOrNqayEPn0MaklSzhjUW+Kz1JKkHCrLdQH5aOikh+rGb17ZiSMMaklSjtii3oL3uvezRS1JyhmDegsW9ugHixbB6tW5LkWSVIQM6i14u9e2mZE338xtIZKkomRQb8G83tmgfv313BYiSSpKBvUWzOs9KDMyd25uC5EkFSWDeguWd+4GffvaopYk5YRB3RTDh9uiliTlhEHdFDvvbItakpQTzQ7qiBgRETPqDcsj4mubLDMuIpbVW+b7LS85B4YPh/nzobo615VIkopMs3smSym9CowGiIhS4B3g3gYW/WtK6bjm7icv7LwzpJR5RGvkyFxXI0kqIq116vtI4I2U0luttL38Mnx45tPr1JKkdtZaQX0acGcj8w6KiBcj4uGI2L2xDUTExIiYFhHTFi9e3EpltZKdd858ep1aktTOWhzUEdEJOB74XQOz/wHskFLaC/gf4L7GtpNSuiGlNCalNKZ///4tLat19e6ded2lLWpJUjtrjRb1scA/UkrvbzojpbQ8pbQyOz4ZKI+Ifq2wz/Y3fLgtaklSu2uNoD6dRk57R8Q2ERHZ8f2z+1vSCvtsfzvvbItaktTuWhTUEVEJfBr4Q71pF0bEhdmvnwdmRcSLwLXAaSml1JJ95oyPaEmScqDZj2cBpJSqgL6bTLu+3vh1wHUt2Ufe8BEtSVIO2DNZU9U+ouV1aklSOzKom6r2ES2vU0uS2pFB3VR9+mQGW9SSpHZkUG8N7/yWJLUzg3pr+Cy1JKmdGdRbY+ed4e23fURLktRuDOqtMXx45hGtf/4z15VIkoqEQb01vPNbktTODOqt4bPUkqR2ZlBvjT59Mm/SskUtSWonBvXW8s5vSVI7Mqi3ls9SS5LakUG9tYYPzzyitWZNriuRJBUBg3prjRiReUTrtddyXYkkqQgY1Ftr1KjM54sv5rYOSVJRMKi31ogR0KmTQS1JahcG9dYqK4PddzeoJUntwqBujr32MqglSe2iLNcFFIKhkx7a6Pt5b5fz/UWLGPNvt/NB197Mu3JCjiqTJHV0tqibYc6AYQDstsiXc0iS2pZB3QyzDWpJUjsxqJthWZfuLOzez6CWJLU5g7qZ5gwYZlBLktqcQd1McwYMY6elC+i0fl2uS5EkdWAGdTPN6T+M8poNDF/ydq5LkSR1YAZ1M3nntySpPRjUzTSv97asLqswqCVJbcqgbqaaklJe7b+9QS1JalMGdQvM6T+MXRfPy7z2UpKkNmBQt8CcAcPos3o5LFyY61IkSR2UQd0CtTeU+YIOSVJbMahb4BWDWpLUxgzqFlhR0ZX5PQca1JKkNmNQt9CcAcMMaklSmzGoW2hO/2Hw2muwenWuS5EkdUAGdQvNHjAMampg5sxclyJJ6oAM6hZ6YdCIzMizz+a2EElSh2RQt9Ci7n1h6FB45plclyJJ6oAM6tZwyCGZoLaHMklSKzOoW8PYsfDee/Dmm7muRJLUwRjUrWHs2Mynp78lSa3MoG4Nu+8OPXsa1JKkVtfioI6IeRHxUkTMiIhpDcyPiLg2Il6PiJkRsU9L95l3Skrg4IPh6adzXYkkqYNprRb14Sml0SmlMQ3MOxYYnh0mAr9spX3ml0MOgdmzYenSXFciSepA2uPU9wnAb1LGc0CviNi2HfbbvmqvU0+dmts6JEkdSmsEdQIejYjpETGxgfmDgfn1vi/ITttIREyMiGkRMW3x4sWtUFY7228/KCvz9LckqVW1RlCPTSntQ+YU90URcegm86OBdT7xwHFK6YaU0piU0pj+/fu3QlntrLIS9t3XG8okSa2qxUGdUlqY/VwE3Avsv8kiC4Dt6n0fAixs6X7z0tix8Pe/w5o1ua5EktRBtCioI6JrRHSvHQeOBmZtstgDwBezd38fCCxLKb3bkv3mrbFjoboa/vGPXFciSeogylq4/kDg3oio3dZvU0qPRMSFACml64HJwHjgdaAKOLeF+8xf9Ts+Oeig3NYiSeoQWhTUKaU3gb0amH59vfEEXNSS/RSMgQNh550zQf2Nb+S6GklSB2DPZK3NF3RIklqRQd3axo6FxYth7txcVyJJ6gAM6tb2qU9lPv/859zWIUnqEAzq1rbLLjBsGDz8cK4rkSR1AAZ1a4uA8ePhsccyj2pJktQCBnVbGD8eqqrgqadyXYkkqcC19DlqAUMnPbTR93nfPxw6d4bJk+Hoo3NUlSSpI7BF3RYqK+GIIzJBLUlSCxjUbWX8+MwjWj6mJUlqAYO6rRx7bObTVrUkqQUM6ray446w664GtSSpRQzqtjR+PDz5JKxaletKJEkFyqBuS+PHw9q19lImSWo2g7otHXIIdOvm6W9JUrMZ1G2pogKOOgoeesi3aUmSmsWgbmvjx8P8+fDyy7muRJJUgAzqtjZ+fObzvvtyW4ckqSDZhWgb+ESXop/6FNxxB3z3u5mXdkiS1ES2qNvDmWfCK6/ACy/kuhJJUoExqNvD5z8P5eVw++25rkSSVGAM6vbQpw9MmAB33gkbNuS6GklSATGo28uZZ8J779n5iSRpqxjU7WDopIcYMbWE5RVd+f2lP851OZKkAmJQt5M1ZZ2YPGIsx7z2LFRV5bocSVKBMKjb0f0jx9Ft7Wp48MFclyJJKhAGdTt6bvs9WNi9n3d/S5KazKBuRylKeGDkYfDII/DBB7kuR5JUAAzqdnbfyHGwfj3cdVeuS5EkFQC7EG1nrwwYxsxtdqbysh/z6fk7kKKEeVdOyHVZkqQ8ZYs6B24ecwI7L13Ap/5pl6KSpM0zqHPgoV0PYVHX3pw37YFclyJJynOe+s6BdaXl/GafCXzjr7ez0wfzP/G2LcDT4ZIkwBZ1zvx29LGsKS3n3Om2qiVJjTOoc2RpZU/uGzmOz836Mz1Xr8h1OZKkPGVQ59AtY46ny/o1nP7in3JdiiQpTxnUOfTKgGE8s8MovviPP1K2YX2uy5Ek5SGDOsdu2fcEBq34IPOyDkmSNmFQ59ifdxrDm70HcdHU/yNSTa7LkSTlGYM6x2pKSvn5IWew2+J5HD/7L7kuR5KUZwzqPPDgbp9i9oBhXPr0HZRvWJfrciRJecSgzgMpSrjq0C+yw0fvcerMKbkuR5KURwzqPPHkjmN4fshIvvLMnXReV53rciRJeaLZQR0R20XEExExJyJejoivNrDMuIhYFhEzssP3W1ZuBxbBVYedzYBVH3Lu9AdzXY0kKU+0pK/v9cDXU0r/iIjuwPSImJJSmr3Jcn9NKR3Xgv0UjWlDdufxnfbjwud+Dx/+DHr3znVJkqQca3aLOqX0bkrpH9nxFcAcYHBrFVasrj70LHquWQU//nGuS5Ek5YFWuUYdEUOBvYG/NTD7oIh4MSIejojdW2N/HdmcATvy+z2OhJ/+FF5+OdflSJJyrMVBHRHdgHuAr6WUlm8y+x/ADimlvYD/Ae7bzHYmRsS0iJi2ePHilpZV0K44/Dzo0QMmToQaO0GRpGLWoqCOiHIyIX1HSukPm85PKS1PKa3Mjk8GyiOiX0PbSindkFIak1Ia079//5aUVfCWVvbk0gO/CM8+y/eO/bcG31ctSSoOLbnrO4CbgDkppZ82ssw22eWIiP2z+1vS3H0Wkz/sfgR/3WE0//7krQxc8UGuy5Ek5UhLWtRjgbOAI+o9fjU+Ii6MiAuzy3wemBURLwLXAqellFILay4OEXz3MxfRqWY9lz92Q66rkSTlSLMfz0opPQ3EFpa5Driuufsodm/33pafj/0Ck/5yK9x3H5x4Yq5LkiS1M3smy3O/2u9EZg8YBhdeCO+9l+tyJEntzKDOc+tLy/jqcd+A5cvhC1+ADRtyXZIkqR0Z1AVgbv8d4Je/hCeegB/8INflSJLakUFdKM4+G849F/7rv+DRR3NdjSSpnRjUheS662D33eGMM+Cdd3JdjSSpHRjUhaSyEn73O1i9Gk45Bap9HaYkdXQGdaHZdVe45RZ49lk46yy7GJWkDq4lr7lUO9q4G9FKvnT4l/iP39/ELfufyA+OnMi8Hx+3meVh3pUT2qFKSVJrM6gL1E37f5ZtV3zA+dPu593u/QBf+S1JHZFBXcB+dMSXGLhyKd958ha4/Qg488xGl7WFLUmFyaAuYClK+PqES+lbtYyDzzkHysvh1FNzXZYkqRV5M1mBW1tWzgUnfQ/Gjs30XHbLLbkuSZLUimxRdwCrKirh4Yfhs5+F886DqipgaK7LkiS1AlvUHUVlJTzwAJxwAlx8MRc+93vwjaKSVPAM6o6koiLTIcrppzPpL7fyo0d/QfmGdbmuSpLUAgZ1R1NeDrfdxi8OPJkzZjzC7Xd9jz5Vy3JdlSSpmbxG3UFs+vgVh53Nq/2HctXD1/DAry/hgs99jzkDdsxNcZKkZrNF3YE9MPIwTv7Cjymt2cA9t3+TU1581OvWklRgDOoO7qVth3P82T/jhUEjuOqRa7n+vivotXp5rsuSJDWRQV0EFnfrw5mn/hf/dfh5HP7G3/nTzRf7TmtJKhAGdZFIUcKv9j+JE7/4U5ZVdIPPfAa+9CX44INclyZJ2gxvJisycwbsyP87+2e8Wvoc/OxncN99cNVVcO65UPLx321t3Te4fY9LUtPYoi5Ca8orMuH8wguw++5w/vnwqU/BtGm5Lk2StAmDukgNnfQQQ29/i6EH/TtfH38JzJ0L++0HJ58Mr76a6/IkSVkGdbGL4J49j4Q33oDLLoNHHoHdd+eKR/6HQcsX5bo6SSp6XqNWRvfucPnl8K//Cj/6EZ+/7n85+aXHeGDkYVy//+e2uLrXnCWpbRjU2tiAAXDNNYxbN5rzn7+P02b+ic/N+jMsehi+/nU47DCIyHWVklQ0PPWtBi3sMYD/PGoiB//LLfxs7Bdg6lQ4/PDMzWfXXgsffZTrEiWpKBjU2qyPuvTgmkO+APPnwy23ZE6Rf/WrMHgwnHMOPP44bNiQ6zIlqcMyqNU0Xbpkgvlvf4Pp0+GMM+Dee+Goo2CHHZj0xM2MfP9N+xKXpFbmNWoBDbx9awvz591wA1xzDTz4INx2G1+afD8XPv8H3uq1DY/scjBM7QMHHLBRJyotqac9bk7zhjhJ+cgWtZqvSxc45RR48EEOuOg3TPrMxbzZZzDnTnsADj4YBg3KtMLvvhs+/DDX1UpSQbJFrVaxtLInd40+hrtGH0P3Nat4acy6TGv7gQfg17/OtKwPOACOOAIOP5yKdWsyPaS1IVvIkjoCg1qtbkVFVzhjQuY69vr18Pzz8PDDMGUKXHkl/OhHzCwtY8a2I5g2ZCTTB+8GSw6Evn1zXbok5R2DWm2rrCxzGvzgg+GHP4Tly+Hpp/n15b/igPmzmPj8Hyiv2QD3/Cdv9BnCi9sO56VtdmafbXZh9sBhVJd3zvVPIEk5ZVCrffXoAePHc8VTmbvDO6+rZtR7r7PvO3PY551XOPitFznp5ScA2BAl/LP3IGYP3BF6vQR77QV77AFDhtjpiqSiYVArp6rLO/P8dnvw/HZ71E0bsGIJe77/OqPefZ2Ri95kn3degW8/9fFK3bvDyJGZYdddYZddYMQI2HFHqGjb696S1N4MarWJltzItah7Xx7v3pfHdz7g4/UnjYWZM2H2bHj55czn5MmZTlhqlZTA9tvDTjvBTjvx5dfW8nbPbXin5wAW9ByYeca7Xkt8S4+kSVI+MKhVGHr1gkMPzQz1LVsGr72WGV59FV5/Hd58E/7wB779wQcbL3vTBZkg3247GDKES16tZmGPfrzfrQ/vd+/L+9368mGX7qTwqUVJ+SNSHvYkNWbMmDRt2rRW254tp+Kwaat9z0v+j8HLFjFk2SKGLHufy0d1hbffznSHOn8+Ne++Rwkb//6vLSnjg669WNy1N4uznx907c2Syp5cfv7h0L9/ZujXL3OXeiueavdxMql4RcT0lNKYhubZolaHtaKiK68MGMYrA4YBcPkmwTfim/cxYOWHDFy5hAErlzJw5VIGrlxC/5Uf0X/Vh2y7Ygmj3nudPlXLKEs18PgNn9xJ9+6ZwK4/9O79yaFXr8zQs+fHQ5n/+0naMv+lUNFaV1rOOz0H8E7PAZtdLlINPatX0nfVMvpXfUiv1SvoW7WM3quX06dqOb2qV9D7w+X0XjiPXqtfomf1SnqsWUVpqtl8AZWVmbvge/aE7t25Y/FaVnaqZGWnLlR16gLpqcwfAt26fTx07Zr5rKzMjNcOlZWZnuJKS7fqGBRaK77Q6i0UWzquHvfcMqilLUhRwkddevBRlx68wXZNWidSDd3Wrs6EdvUqeqzJfPasXkm3tVV0X1PFpQdsk7nGvnw5LF9Op3ffZodV79J1XTWVa1fD7MdhzZqtK7Zz549Du/5n586Z8doh+/3bzy9kTWk5a8o6saasHK77Z2ZeRUVmqD/eqdMnx+tPKy/PDD46J7WqFgV1RBwDXAOUAr9KKV25yfwK4DfAvsAS4NSU0ryW7FMqBClKWFHRNdNLW8+Gl7l0k1bJyQ21Wtatg5UrYdWqzGft+KZDVdXHw6pVsHp1Zqiq+nh88eKPx6urobqas1asovO6tR9fq3+iFX748vJMgNd+1o7Xn97YUFb2yfGysrrxbz39FutLSusGfjzr4/mlpZsfLy395Hj9obHppaWZJwo2N72hT/9gUStp9s1kEVEKvAZ8GlgA/B04PaU0u94y/wqMSildGBGnAZ9NKZ26pW17M5m05dOPrbUOKVFWs4GK9Wt5+buHZ4J8zRpYs4bjfvIYndavo2LDOso3rOM3Z46GtWsz89eu/Xh8zZrMHxX1p9V+r/f5yAtvU7ZhPeU1GyirWU/Zhg2U16ynrGYDZTUbGNmvS2b5detg/XoWfbiK0poNlNdsoDS7TGnNhsw9A/kuIhPYDYV4vWFR1XpqIqiJEmoiSNnP2mk7Dei+8Tr1t7vpUH/epuObmTf55fdJ2X0n4P+NHrzROve88E52XqauU/ff/uP5tdva3PfGpmWH6558o27bEKSARGRqAr7+mV0bXbduqD3mm45vadktzWtsu/vtB7vt1oq/Lm1zM9n+wOsppTezO7kLOAGYXW+ZE4DLs+O/B66LiEj5eKu5VKwiWF9axvrSsszd7PXM2mbexsse37Jrkxdu4Q+HTf/Q2L+x5VOiNNVQtuHjkC9NNfzj24fDhg2ZPubXr/94vKHPTYfs9C/f+jylqYbSmg2UpBquOWWvjZb73j0vUpLdX0lK/Mf4XT+eX1OTGWq/p/Tx9Prz631/bOo8Irut0lRDkIiUKEk1lNbUsNOobTdeJ6XMUP97/e029r2hz+z6Oy1dTklKlKQEJJj+3kbL7b+0ikgf18WSlz+uo349m/ve2LSUuHBDDQGZ7dNAPDzdrF+3tvXzn7dqUG9OS1rUnweOSSmdn/1+FnBASuniesvMyi6zIPv9jewyHzSwvYnAxOzXEcCrzSqs4+gHfOI4qdV5nNuHx7n9eKzbR2sf5x1SSv0bmtGSFnVDF2A2Tf2mLJOZmNINQAPPvxSniJjW2GkQtR6Pc/vwOLcfj3X7aM/j3JIumBbARrfADgEWNrZMRJSRua1maQv2KUlSUWlJUP8dGB4RwyKiE3Aa8MAmyzwAnJ0d/zzwZ69PS5LUdM0+9Z1SWh8RFwN/IvN41s0ppZcj4j+BaSmlB4CbgNsi4nUyLenTWqPoIuFlgPbhcW4fHuf247FuH+12nPOyr29JkpTha4IkScpjBrUkSXnMoM6xiDgmIl6NiNcjYlID8y+NiNkRMTMiHo+IHXJRZ6Hb0nGut9znIyJFhI+3NENTjnNEnJL9nX45In7b3jV2BE34d2P7iHgiIl7I/tsxPhd1FrqIuDkiFmX7BGlofkTEtdn/DjMjYp82KSSl5JCjgcxNeG8AOwKdgBeBkZssczhQmR3/F+DuXNddaENTjnN2ue7AU8BzwJhc111oQxN/n4cDLwC9s98H5LruQhuaeJxvAP4lOz4SmJfrugtxAA4F9gFmNTJ/PPAwmT5DDgT+1hZ12KLOrbpuWFNKa4HabljrpJSeSClVZb8+R+Z5dW2dLR7nrB8CVwHV7VlcB9KU43wB8IuU0ocAKaVF7VxjR9CU45yAHtnxnnyyjws1QUrpKTbf98cJwG9SxnNAr4jYtrXrMKhzazAwv973BdlpjfkSmb/etHW2eJwjYm9gu5TSH9uzsA6mKb/PuwC7RMQzEfFc9g182jpNOc6XA2dGxAJgMvBv7VNa0dnaf8ObxfdR51aTu1iNiDOBMcBhbVpRx7TZ4xwRJcDPgHPaq6AOqim/z2VkTn+PI3N26K8RsUdK6aM2rq0jacpxPh24NaX03xFxEJn+LPZIqRBeO1ZQmvxveEvYos6tpnTDSkQcBXwXOD6ltKadautItnScuwN7AE9GxDwy15oe8IayrdbUboXvTymtSyn9k8zLd4a3U30dRVOO85eA/wNIKU0FOpN5iYRaV5P+DW8pgzq3ttgNa/aU7P9HJqS9ntc8mz3OKaVlKaV+KaWhKaWhZO4FOD6l1HovRS8OTelW+D4yN0gSEf3InAp/s12rLHxNOc5vA0cCRMRuZIJ6cbtWWRweAL6Yvfv7QGBZSund1t6Jp75zKDWtG9afAN2A30XmxeVvp5SOz1nRBaiJx1kt1MTj/Cfg6IiYDWwAvplSWpK7qgtPE4/z14EbI+ISMqdiz0nZ25TVdBFxJ5nLNP2y1/svA8oBUkrXk7n+Px54HagCzm2TOvxvJ0lS/vLUtyRJecygliQpjxnUkiTlMYNakqQ8ZlBLkpTHDGpJkvKYQS1JUh77/wF6ZO/gh9NgFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#check sin(theta/2)\n",
    "x = sampling(func, range_, 1000)\n",
    "ax.cla()\n",
    "ax.plot(x0, px0, 'r-', label='Theory')\n",
    "ax.hist(x, bins=100, range=(0.1, 1), density=True, label='Monte-Carlo Sampling Result')\n",
    "ax.legend()\n",
    "ax.set_title(r'$\\sin\\theta/2$')\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAHiCAYAAAA9Am/ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAXfElEQVR4nO3df7Dsd13f8dfbJNIfaIPmQFNIesEJTtHRQG8zdBgtLbYG6BBtxSbTgUCxV1todfSPRuyIdcaZ2IpOrS1MKBlCByLUQIkmWlOkMnYKeoMxBgMaaJQrmeRKLNAJpZP47h9nL55cTnIOd8+55312H4+ZnbP72e/ufnbZ8Dyf7+753uruAAAzfNlBTwAA+FPCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIsywpqrq2VX1P6rqoar69aq6+KDnBAgzrKWqelqSW5P8eJKvTvLxJP/yQCcFJBFmWFevT/Km7r65uz+X5GeT/LUDnhOQ5NyDngBwdlXVVya5Iskztwx/WZL/ezAzArYSZlg/L0hyXpI7q+rU2BOSvOfAZgR8gV3ZsH6OJLm5u88/dUryviS/dLDTAhJhhnX0hCQPnbpQVU9PcjTJzQc2I+ALhBnWz28k+RtV9Zeq6qIkb0/yQ9394AHPC4jPmGEd/UqSn0/yu0k+leTHu/tNBzsl4JTq7oOeAwCwYFc2AAwizAAwiDADwCDCDACDCDMADDLiz6UuuOCCPnLkyEFPAwDOmttvv/2Punvj9PERYT5y5EiOHz9+0NMAgLOmqn5/u3G7sgFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAY596AnsCqOXHPLoy7fe+2LD2gmABxmVswAMIgwA8AgwgwAgwgzAAwizAAwyI5hrqqLqup9VXV3VX24qr53Mf5VVXVbVf3e4ueTFuNVVT9dVfdU1Z1V9Zz9fhIAsCp2s2J+OMkPdPdfSfLcJK+uqmcluSbJe7v7kiTvXVxOkhcmuWRxOpbkDXs+awBYUTuGubvv6+4PLc5/NsndSZ6a5IokNyw2uyHJty3OX5Hkrb3pA0nOr6oL93zmALCCvqTPmKvqSJJnJ/lgkqd0933JZryTPHmx2VOTfGLLzU4sxgCAHew6zFX1xCQ3Jfm+7v7M4226zVhvc3/Hqup4VR0/efLkbqcBACttV2GuqvOyGeW3dfe7FsP3n9pFvfj5wGL8RJKLttz8aUk+efp9dvd13X20u49ubGyc6fwBYKXs5lvZleTNSe7u7p/cctXNSa5enL86yXu2jL988e3s5yb59Kld3gDA49vNP2LxvCQvS/LbVXXHYuy1Sa5N8s6qelWSP0jy0sV1tyZ5UZJ7kjyU5JV7OmMAWGE7hrm7fy3bf26cJC/YZvtO8uol5wUAa8mRvwBgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQXYMc1VdX1UPVNVdW8beUVV3LE73VtUdi/EjVfW5Lde9cT8nDwCr5txdbPOWJD+T5K2nBrr7H5w6X1WvT/LpLdt/rLsv3asJTnXkmlsOegoArKAdw9zd76+qI9tdV1WV5DuT/K29nRYArKdlP2P+piT3d/fvbRl7elX9ZlX9alV902PdsKqOVdXxqjp+8uTJJacBAKth2TBfleTGLZfvS3Jxdz87yfcneXtVfeV2N+zu67r7aHcf3djYWHIaALAazjjMVXVukr+X5B2nxrr78939qcX525N8LMkzl50kAKyLZVbM35LkI9194tRAVW1U1TmL889IckmSjy83RQBYH7v5c6kbk/zPJF9bVSeq6lWLq67Mo3djJ8k3J7mzqn4ryc8l+Z7ufnAvJwwAq2w338q+6jHGX7HN2E1Jblp+WgCwnhz5CwAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAbZMcxVdX1VPVBVd20Z+5Gq+sOqumNxetGW636wqu6pqo9W1bfu18QBYBXtZsX8liSXbzP+U9196eJ0a5JU1bOSXJnk6xa3+Q9Vdc5eTRYAVt2OYe7u9yd5cJf3d0WSn+3uz3f3/0pyT5LLlpgfAKyVZT5jfk1V3bnY1f2kxdhTk3xiyzYnFmMAwC6caZjfkORrklya5L4kr1+M1zbb9nZ3UFXHqup4VR0/efLkGU4DAFbLGYW5u+/v7ke6+0+SvCl/urv6RJKLtmz6tCSffIz7uK67j3b30Y2NjTOZBgCsnDMKc1VduOXityc59Y3tm5NcWVVPqKqnJ7kkya8vN0UAWB/n7rRBVd2Y5PlJLqiqE0lel+T5VXVpNndT35vku5Okuz9cVe9M8jtJHk7y6u5+ZH+mDgCrZ8cwd/dV2wy/+XG2/7EkP7bMpABgXTnyFwAMIswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgwgwAgwgzAAwizAAwyI5hrqrrq+qBqrpry9i/qaqPVNWdVfXuqjp/MX6kqj5XVXcsTm/cz8kDwKrZzYr5LUkuP23stiRf393fkOR3k/zglus+1t2XLk7fszfTBID1sGOYu/v9SR48beyXu/vhxcUPJHnaPswNANbOXnzG/I+S/OKWy0+vqt+sql+tqm96rBtV1bGqOl5Vx0+ePLkH0wCAw2+pMFfVDyV5OMnbFkP3Jbm4u5+d5PuTvL2qvnK723b3dd19tLuPbmxsLDMNAFgZZxzmqro6yd9N8g+7u5Okuz/f3Z9anL89yceSPHMvJgoA6+CMwlxVlyf5F0le0t0PbRnfqKpzFuefkeSSJB/fi4kCwDo4d6cNqurGJM9PckFVnUjyumx+C/sJSW6rqiT5wOIb2N+c5Eer6uEkjyT5nu5+cNs7XnFHrrnlUZfvvfbFBzQTAA6THcPc3VdtM/zmx9j2piQ3LTspAFhXjvwFAIMIMwAMIswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgwgwAg+wqzFV1fVU9UFV3bRn7qqq6rap+b/HzSYvxqqqfrqp7qurOqnrOfk0eAFbNblfMb0ly+Wlj1yR5b3dfkuS9i8tJ8sIklyxOx5K8YflpAsB62FWYu/v9SR48bfiKJDcszt+Q5Nu2jL+1N30gyflVdeFeTBYAVt0ynzE/pbvvS5LFzycvxp+a5BNbtjuxGAMAdrAfX/6qbcb6izaqOlZVx6vq+MmTJ/dhGgBw+CwT5vtP7aJe/HxgMX4iyUVbtntakk+efuPuvq67j3b30Y2NjSWmAQCrY5kw35zk6sX5q5O8Z8v4yxffzn5ukk+f2uUNADy+c3ezUVXdmOT5SS6oqhNJXpfk2iTvrKpXJfmDJC9dbH5rkhcluSfJQ0leucdzBoCVtaswd/dVj3HVC7bZtpO8eplJAcC6cuQvABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBdnWsbJZ35JpbHnX53mtffEAzAWAyK2YAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYJBzz/SGVfW1Sd6xZegZSX44yflJ/nGSk4vx13b3rWc8QwBYI2cc5u7+aJJLk6Sqzknyh0neneSVSX6qu39iT2YIAGtkr3ZlvyDJx7r79/fo/gBgLe1VmK9McuOWy6+pqjur6vqqetJ2N6iqY1V1vKqOnzx5crtNAGDtLB3mqvryJC9J8p8XQ29I8jXZ3M19X5LXb3e77r6uu49299GNjY1lpwEAK2EvVswvTPKh7r4/Sbr7/u5+pLv/JMmbkly2B48BAGthL8J8Vbbsxq6qC7dc9+1J7tqDxwCAtXDG38pOkqr6c0n+dpLv3jL8r6vq0iSd5N7TrgMAHsdSYe7uh5J89WljL1tqRgCwxhz5CwAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBljrACGfuyDW3POryvde++IBmAsAkVswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDnHvQE2HTkmlsedfnea198QDMB4CBZMQPAIEuvmKvq3iSfTfJIkoe7+2hVfVWSdyQ5kuTeJN/Z3X+87GMBwKrbqxXz3+zuS7v76OLyNUne292XJHnv4jIAsIP92pV9RZIbFudvSPJt+/Q4ALBS9iLMneSXq+r2qjq2GHtKd9+XJIufT96DxwGAlbcX38p+Xnd/sqqenOS2qvrIbm60iPixJLn44ov3YBoAcPgtvWLu7k8ufj6Q5N1JLktyf1VdmCSLnw9sc7vruvtodx/d2NhYdhoAsBKWCnNV/fmq+opT55P8nSR3Jbk5ydWLza5O8p5lHgcA1sWyu7KfkuTdVXXqvt7e3b9UVb+R5J1V9aokf5DkpUs+DgCshaXC3N0fT/KN24x/KskLlrnvdedIYADryZG/AGAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAYRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQYQZAAY596AnwO4cueaWR12+99oXH9BMANhPVswAMIgwA8AgwgwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIMIMwAMIswAMIgwA8AgZxzmqrqoqt5XVXdX1Yer6nsX4z9SVX9YVXcsTi/au+kCwGpb5t9jfjjJD3T3h6rqK5LcXlW3La77qe7+ieWnx26d/u81J/7NZoDD6IzD3N33Jblvcf6zVXV3kqfu1cQAYB0ts2L+gqo6kuTZST6Y5HlJXlNVL09yPJur6j/e5jbHkhxLkosvvngvprFWtlshA3D4Lf3lr6p6YpKbknxfd38myRuSfE2SS7O5on79drfr7uu6+2h3H93Y2Fh2GgCwEpYKc1Wdl80ov62735Uk3X1/dz/S3X+S5E1JLlt+mgCwHpb5VnYleXOSu7v7J7eMX7hls29PcteZTw8A1ssynzE/L8nLkvx2Vd2xGHttkquq6tIkneTeJN+91AwBYI0s863sX0tS21x165lPBwDWmyN/AcAgwgwAgwgzAAwizAAwyJ4c+YuZTj86mGNnA8xnxQwAgwgzAAwizAAwiDADwCDCDACDCDMADCLMADCIMAPAIMIMAIM48tcacSQwgPmsmAFgEGEGgEGEGQAGEWYAGESYAWAQYQaAQfy5FF/wpf45lT+/Ath7VswAMIgVM7t2+goZgL1nxQwAgwgzAAwizAAwiM+Y15jPjAHmsWIGgEGsmBnD30UDWDEDwChWzBwYn3EDfDErZgAYxIqZx7TXK1orZICdWTEDwCBWzOwZ36oGWJ4VMwAMYsXMobHTityKHVgFVswAMIgVM/vGt7ABvnRWzAAwyL6tmKvq8iT/Nsk5Sf5jd1+7X4/Fatppxb0fK/Kd7vNsf27tc/OdeY1YNfuyYq6qc5L8+yQvTPKsJFdV1bP247EAYJXs14r5siT3dPfHk6SqfjbJFUl+Z58eD77I2VhRf6mrs2XndDYef6/v86BXsAcxP39BwDL26zPmpyb5xJbLJxZjAMDjqO7e+zutemmSb+3u71pcflmSy7r7n23Z5liSY4uLX5vko3s8jQuS/NEe3+dhsu7PP/EarPvzT7wGiddg8vP/y929cfrgfu3KPpHkoi2Xn5bkk1s36O7rkly3T4+fqjre3Uf36/6nW/fnn3gN1v35J16DxGtwGJ//fu3K/o0kl1TV06vqy5NcmeTmfXosAFgZ+7Ji7u6Hq+o1Sf5rNv9c6vru/vB+PBYArJJ9+zvm7r41ya37df+7sG+7yQ+JdX/+iddg3Z9/4jVIvAaH7vnvy5e/AIAz45CcADDIoQ5zVV1eVR+tqnuq6pptrn9CVb1jcf0Hq+rI2Z/l/trFa/CKqjpZVXcsTt91EPPcL1V1fVU9UFV3Pcb1VVU/vXh97qyq55ztOe63XbwGz6+qT295D/zw2Z7jfqqqi6rqfVV1d1V9uKq+d5ttVvZ9sMvnv+rvgT9TVb9eVb+1eA3+1TbbHJ4edPehPGXzS2UfS/KMJF+e5LeSPOu0bf5pkjcuzl+Z5B0HPe8DeA1ekeRnDnqu+/gafHOS5yS56zGuf1GSX0xSSZ6b5IMHPecDeA2en+QXDnqe+/j8L0zynMX5r0jyu9v8d7Cy74NdPv9Vfw9Ukicuzp+X5INJnnvaNoemB4d5xfyFw3529/9Lcuqwn1tdkeSGxfmfS/KCqqqzOMf9tpvXYKV19/uTPPg4m1yR5K296QNJzq+qC8/O7M6OXbwGK6277+vuDy3OfzbJ3fniIw2u7Ptgl89/pS3+d/0/i4vnLU6nf4Hq0PTgMId5N4f9/MI23f1wkk8n+eqzMruzY7eHPv37i913P1dVF21z/SpzeNhNf32xm+8Xq+rrDnoy+2Wxe/LZ2VwxbbUW74PHef7Jir8HquqcqrojyQNJbuvux3wPTO/BYQ7zdr/pnP4b0m62Ocx28/x+PsmR7v6GJP8tf/ob47pY9ffAbnwom4f++8Yk/y7Jfzng+eyLqnpikpuSfF93f+b0q7e5yUq9D3Z4/iv/HujuR7r70mweafKyqvr60zY5NO+BwxzmHQ/7uXWbqjo3yV/Iau3y282hTz/V3Z9fXHxTkr96luY2xW7eJyutuz9zajdfbx5f4LyquuCAp7Wnquq8bEbpbd39rm02Wen3wU7Pfx3eA6d09/9O8t+TXH7aVYemB4c5zLs57OfNSa5enP+OJL/Si0/+V8SOr8Fpn6O9JJufP62Tm5O8fPGt3Ocm+XR333fQkzqbquovnvosraouy+Z/95862FntncVze3OSu7v7Jx9js5V9H+zm+a/Be2Cjqs5fnP+zSb4lyUdO2+zQ9GDfjvy13/oxDvtZVT+a5Hh335zNN+t/qqp7svmb0ZUHN+O9t8vX4J9X1UuSPJzN1+AVBzbhfVBVN2bzG6cXVNWJJK/L5hc/0t1vzObR516U5J4kDyV55cHMdP/s4jX4jiT/pKoeTvK5JFdO/T+kM/S8JC9L8tuLzxiT5LVJLk7W4n2wm+e/6u+BC5PcUFXnZPOXjnd29y8c1h448hcADHKYd2UDwMoRZgAYRJgBYBBhBoBBhBkABhFmABhEmAFgEGEGgEH+PyQZmlcnoRujAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#watch theta\n",
    "theta = np.arcsin(x) * 2\n",
    "ax.cla()\n",
    "ax.hist(theta, bins=100, range=(0, np.pi), density=False)\n",
    "ax.set_title(r'$\\theta$')\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 反函数直接抽样法\n",
    "\n",
    "由于$\\theta \\in (0, \\pi)$，显然$\\theta/x\\in (0, \\pi/2)$，于是$\\sin(\\theta/2)\\in(0, 1)$\n",
    "\n",
    "也即上面的$x\\in(0,1)$\n",
    "\n",
    "由于做了0.1的截断，所以x的范围设定在$(0.1, 1)$\n",
    "\n",
    "十分显然地，有：\n",
    "\n",
    "$$F(x) = \\left\\{ \\begin{aligned}0, x<0.1;\\\\\n",
    "A(-\\frac{1}{2}x^{-2}+\\frac{1}{2}(-0.1)^{-2}), 0.1\\le x \\le 1;\\\\\n",
    "1\\end{aligned}\\right.$$\n",
    "\n",
    "注意到两头的并不需要我们管，因为我们对cdf的抽样在0到1，它们可以看做分别对应了$0-\\epsilon$和$1+\\epsilon$\n",
    "\n",
    "所以我们关注中间这一坨就行了，取$x = 1$的时候，对应的就是原函数在$(0.1, 1)$的定积分结果（不考虑归一化），很容易算出来：\n",
    "\n",
    "$A^{-1} = -\\frac{1}{2}(1^{-2}-(0.1)^{-2})\\Rightarrow A = 2/99$\n",
    "\n",
    "所以有$$F(x) = \\frac{1}{99}(100 - x^{-2}), 0.1\\le x \\le 1$$\n",
    "\n",
    "显然，反函数为：\n",
    "\n",
    "$$ x(F) = (100 - 99F)^{-1/2} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#inverse sampling\n",
    "def inverse_sampling(func, n=1000):\n",
    "    r'''\n",
    "    Inversed function sampling.\n",
    "    Args: 2\n",
    "        func: a function, the inverse function. you need to derive its analytical form in advance.\n",
    "        n: an integer, how many variables you want to simulate.\n",
    "    Return: 1\n",
    "        res: an array with n variables respect to the inversed function of func (pdf.).\n",
    "    '''\n",
    "    y = np.random.uniform(0, 1, (n,))\n",
    "    x = func(y)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAHjCAYAAADsYlzyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU9b3/8dcn+0JkDQgECVgRQSHaiAug4A6idHGtu7WIV39WutL2XrVWW6xY116pu9dSxa2KotYN6gbFCBgBERCjRlTCGkIWSPL9/TGTGMKETDKTOTM57+fjcR5z5pwzZz450bz5fs8532POOURERCQ+JXldgIiIiLRMQS0iIhLHFNQiIiJxTEEtIiISxxTUIiIicUxBLSIiEscU1CIiInFMQS0iIhLHFNQinZyZrTCzcV7XISLto6AW6eScc8OdcwvC3d7MDjWzd8ys0swWm9l+Tdb1M7NSM0s3swfM7DMz225mS81sQof8ACI+p6AWkUZmlge8CNwM9ATWAf/dZJOJwMtACvAFcCzQFfgf4Akzy49huSK+oKAW6STM7Ndm9mWwhfuxmR0fXF5iZic0mf+FmRWb2TYzm2NmGU12cytwn3NurnOuCngcOLzJ+onAi865Hc65651zJc65eufcC8CnwHdj89OK+IeCWqQTMLMDgauAw51zOcDJQEkLm58FnAIMAkYAFwf3sQ8wGbi/ybZJQHVwfSpwDPBqiO/vAwwBVkT8w4jIblK8LkBEoqIOSAeGmVmZc65kL9ve6ZxbD2BmzwMFweXHA6lAsZk1bJsOPBecPwb4wDm3venOggE+G3jEObcqCj+LiDShFrVIJ+CcWwtcA1wPbDCzx82sXwubf91kvhLoEpzPB+Y657o1TMB8AuekIdjt3XRHZpYEPArsJNCiF5EoU1CLdBLOuX8458YAAwFH4IKwtkgnENwAmNkgoBCYG1w0EZjXZL0BDwB9gB8653a1v3oRaYmCWqQTMLMDzew4M0sncE65ikB3eFu8BxwbvAVrAPAP4HfOuc3B0E5v1rV9D3AQcFrwwjMR6QAKapHOIR2YAWwk0LXdG/htG/fxBvA8sBp4G3jUOXdfcN2pNOn2NrOBwOUEzm9/bWYVwem8iH4KEdmDOee8rkFE4pyZvQjc7Zx7sdWNRSSq1KIWkXAsIHBhmYjEmFrUIiIicUwtahERkTimoBYREYljcTkyWa9evVx+fr7XZYiIiMTE+++/v9E5lxtqXVwGdX5+PkVFRV6XISIiEhNm9llL69T1LSIiEscU1CIiInFMQS0iIhLH4vIctYgkpl27dlFaWkp1dbXXpYjEpYyMDPLy8khNTQ37MwpqEYma0tJScnJyyM/Pp8kzrUUEcM6xadMmSktLGTRoUNifU9e3iERNdXU1PXv2VEiLhGBm9OzZs809TgpqEYkqhbRIy9rz/4eCWkREJI7pHLWIdBqbNm3i+OOPB+Drr78mOTmZ3NxcSkpK6NevHytXrvS4QpG2U4taRDqNnj17smzZMpYtW8bUqVOZNm1a4/ukpOj/uautrY36PkWaU4taRDrGNdfAsmXR3WdBAdx+e7s+WldXx09+8hPeffdd+vfvz3PPPUdmZiaffPIJV155JWVlZWRlZXHfffcxdOhQPvvsMy699FLKysrIzc3loYceYr/99uPiiy+mR48eLF26lIKCAl544QXeffddcnNzqa+vZ8iQISxatIhevXpF92cX31KLWkR8Yc2aNVx55ZWsWLGCbt268fTTTwMwZcoU7rrrLt5//31mzpzJf/3XfwFw1VVXceGFF1JcXMx5553H1Vdf3biv1atX89prr3Hbbbdx/vnnM3v2bABee+01Ro4cqZCWqFKLWkQ6Rjtbvh1l0KBBFBQUAPDd736XkpISKioqePfddznzzDMbt6upqQFg4cKFPPPMMwBccMEF/OpXv2rc5swzzyQ5ORmASy+9lMmTJ3PNNdfw4IMPcskll8TqRxKfUFCLiC+kp6c3zicnJ1NVVUV9fT3dunVjWRhd9E1vq8nOzm6cHzBgAH369OGNN97gP//5T2PrWiRa1PUtIr61zz77MGjQIJ588kkgMHLUBx98AMDRRx/N448/DsDs2bMZM2ZMi/u57LLLOP/88znrrLMaW9oi0aKgFhFfmz17Ng888AAjR45k+PDhPPfccwDceeedPPTQQ4wYMYJHH32UO+64o8V9nH766VRUVKjbWzqEOee8rmEPhYWFrqioyOsyRKSNPvroIw466CCvy4i5oqIipk2bxltvveV1KZIAQv1/YmbvO+cKQ23f+c9Rr18PKSnQu7fXlYhIJzRjxgzuuecenZuWDtP5u75HjYLf/tbrKkSkk5o+fTqfffbZXs9hi0Si8wd1djbs2OF1FSIiIu3SalCb2QAzm29mH5nZCjP7aXB5DzN71czWBF+7t/D5i4LbrDGzi6L9A7RKQS0iIgksnBZ1LfBz59xBwJHAlWY2DJgOvO6cOwB4Pfh+N2bWA7gOOAIYBVzXUqB3GAW1iIgksFaD2jn3lXNuSXB+O/AR0B+YDDwS3OwR4HshPn4y8KpzbrNzbgvwKnBKNAoPm4JaREQSWJuu+jazfOBQ4D9AH+fcVxAIczMLdVl1f+CLJu9Lg8tC7XsKMAVgv/32a0tZe5edDV9+Gb39iUjY8qfPi+r+Smac2uo2Xbp0oaKiIqrfG03jxo1j5syZFBbufifOW2+9xdSpU0lNTWXhwoVkZma2ab8LFiwgLS2No48+GoBZs2aRlZXFhRdeGFG99fX1XHPNNbzxxhuYGRkZGTzxxBMMGjQoov3uTX5+PkVFRfTq1Yujjz6ad999N+J9LliwgMmTJzN48GCqqqqYNGkSM2fOjEK137r44ouZNGkSZ5xxBrfffjtTpkwhKysr4v2GfTGZmXUBngaucc6Vh/uxEMtC3rjtnLvXOVfonCvMzc0Nt6zWqUUtIlFUV1fXIfudPXs2v/jFL1i2bFmbQxoCQdQ00KZOnRpxSAPMmTOH9evXU1xczIcffsg///lPunXrFvF+wxWNkG4wduxYli5dytKlS3nhhRd45513orbv5m6//XYqKyujsq+wgtrMUgmE9Gzn3DPBxd+YWd/g+r7AhhAfLQUGNHmfB6xvf7ntoKAW8aUFCxYwbtw4zjjjDIYOHcp5552Hc46XXnqJs846a7ftTjvtNABeeeUVjjrqKA477DDOPPPMxpZ5fn4+N9xwA2PGjOHJJ5/kzjvvZNiwYYwYMYJzzjkHgB07dnDppZdy+OGHc+ihhzaOcFZVVcU555zDiBEjOPvss6mqqtqj1vvvv58nnniCG264gfPOO48FCxYwadKkxvVXXXUVDz/8cGMt1113HYcddhiHHHIIq1atoqSkhFmzZnHbbbdRUFDAW2+9xfXXX9/YYly2bBlHHnkkI0aM4Pvf/z5btmwBAq37X//614waNYohQ4aEHLDlq6++om/fvo3P887Ly6N798ClRldccQWFhYUMHz6c6667rvEz+fn5/Pa3v+Woo46isLCQJUuWcPLJJ7P//vsza9asxuN+zDHH8P3vf59hw4YxdepU6uvr9/j+Ll267PX3CfDiiy8ydOhQxowZw9VXX73bsQslMzOTgoICvgz2trb0u1uxYgWjRo2ioKCAESNGsGbNGkpKSjj44IMb9zVz5kyuv/763fZ/5513sn79esaPH8/48eP3Wks4wrnq24AHgI+cc39psmou0HAV90XAcyE+/i/gJDPrHryI7KTgsthRUIv41tKlS7n99ttZuXIl69at45133uHEE09k0aJF7Aj+XZgzZw5nn302Gzdu5MYbb+S1115jyZIlFBYW8pe/fPsnLyMjg7fffptzzjmHGTNmsHTpUoqLixuD56abbuK4447jvffeY/78+fzyl79kx44d3HPPPWRlZVFcXMzvfvc73n///T3qvOyyyzj99NO55ZZbwho4pVevXixZsoQrrriCmTNnkp+fz9SpU5k2bRrLli1j7Nixu21/4YUXcvPNN1NcXMwhhxzC73//+8Z1tbW1LF68mNtvv3235Q3OOussnn/+eQoKCvj5z3/O0qVLG9fddNNNFBUVUVxczL///W+Ki4sb1w0YMICFCxcyduxYLr74Yp566ikWLVrEtdde27jN4sWLufXWW/nwww/55JNPGp9W1pJQv8/q6mouv/xyXnrpJd5++23KyspaPX5btmxhzZo1HHPMMY0/R6jf3axZs/jpT3/KsmXLKCoqIi8vr9V9A1x99dX069eP+fPnM3/+/LA+szfhtKhHAxcAx5nZsuA0EZgBnGhma4ATg+8xs0Izux/AObcZ+APwXnC6IbgsdrKzobISQvxLTUQ6t1GjRpGXl0dSUhIFBQWUlJSQkpLCKaecwvPPP09tbS3z5s1j8uTJLFq0iJUrVzJ69GgKCgp45JFH+Oyzzxr3dfbZZzfOjxgxgvPOO4+///3vpKQELvV55ZVXmDFjBgUFBYwbN47q6mo+//xz3nzzTc4///zGz40YMSLin+sHP/gB8O3jOvdm27ZtbN26lWOPPRaAiy66iDfffDPsfeXl5fHxxx/zpz/9iaSkJI4//nhef/11AJ544gkOO+wwDj30UFasWMHKlSsbP3f66acDcMghh3DEEUeQk5NDbm4uGRkZbN26FQj8fgYPHkxycjLnnnsub7/99l5/llC/z1WrVjF48ODGc+bnnntui59/6623GDFiBPvuuy+TJk1i3333BVr+3R111FH88Y9/5Oabb+azzz5r1ymJaGj1YjLn3NuEPtcMcHyI7YuAy5q8fxB4sL0FRiw7G5yDqqrAvIj4RvNHW9bW1gKB0P3rX/9Kjx49OPzww8nJycE5x4knnshjjz0Wcl9NH205b9483nzzTebOncsf/vAHVqxYgXOOp59+mgMPPHCPzzZ9RGY4UlJSdusGrq6uDvlzNf2Z2iucfaWnpzNhwgQmTJhAnz59ePbZZxk8eDAzZ87kvffeo3v37lx88cW71dmw36SkpN1+D0lJSY3f0/y4tHacQv0+2/K8irFjx/LCCy+wevVqxowZw/e//30KCgpa/N0ddNBBHHHEEcybN4+TTz6Z+++/nyFDhuz1d9MR/DEyGaj7W0QajRs3jiVLlnDfffc1tpSPPPJI3nnnHdauXQtAZWUlq1ev3uOz9fX1fPHFF4wfP54///nPbN26lYqKCk4++WTuuuuuxuBo6CI+5phjGruzly9fvlv3cEsGDhzIypUrqampYdu2bY0t2L3Jyclh+/bteyzv2rUr3bt3bzz//Oijjza2rsOxZMkS1q8PXFpUX19PcXExAwcOpLy8nOzsbLp27co333zDSy+9FPY+GyxevJhPP/2U+vp65syZ065hWIcOHcq6desaewPmzJnT6meGDBnCb37zG26++WaAFn9369atY/DgwVx99dWcfvrpFBcX06dPHzZs2MCmTZuoqanhhRdeCPkdLf0+2qPzP5RDQS3imXBup/JCcnIykyZN4uGHH+aRRwLDQeTm5vLwww9z7rnnUlNTA8CNN97IkCFDdvtsXV0d559/Ptu2bcM5x7Rp0+jWrRv/8z//wzXXXMOIESNwzpGfn88LL7zAFVdcwSWXXMKIESMoKChg1KhRrdY3YMAAzjrrLEaMGMEBBxzAoYce2upnTjvtNM444wyee+457rrrrt3WPfLII0ydOpXKykoGDx7MQw89FO6hYsOGDfzkJz9pPCajRo3iqquuIiMjg0MPPZThw4czePBgRo8eHfY+Gxx11FFMnz6dDz/8sPHCsrbKzMzkf//3fznllFPo1atXWMcXAlfFz5w5k08//bTF392cOXP4+9//TmpqKvvuuy/XXnstqampXHvttRxxxBEMGjSIoUOHhtz/lClTmDBhAn379o34PHXnf8zlk0/CWWfBhx9Ckyv1RCT6/PqYS2m7BQsWMHPmzBZbpG1RUVFBly5dcM5x5ZVXcsABBzBt2rQoVNkx2vqYS3V9i4hIQrvvvvsoKChg+PDhbNu2jcsvv9zrkqJKXd8iIhJz48aNY9y4cVHZ17Rp0+K6BR0ptahFJKri8XSaSLxoz/8fCmoRiZqMjAw2bdqksBYJwTnHpk2byMjIaNPn1PUtIlGTl5dHaWlpWKNDifhRRkZG2COcNVBQi0jUpKamduhTlUT8SF3fIiIicazzB3V6OiQlKahFRCQhdf6gNtMTtEREJGF1/qAGBbWIiCQsBbWIiEgc80VQf1RexyvvrSN/+jzyp8/zuhwREZGw+SKoK1MzyNzZ8c8MFRERiTbfBHXWLgW1iIgkHl8EdVWaglpERBKTL4K6MjWdzF01XpchIiLSZj4JarWoRUQkMfkiqKtSM8hUUIuISALyRVAHWtQ1oEfviYhIgvFFUFelppPs6kmv2+V1KSIiIm3ii6CuTA08pFvd3yIikmh8FdS6oExERBKNL4K6KjUdgMydukVLREQSiy+CujJNLWoREUlM/ghqdX2LiEiC8kVQV+liMhERSVC+COrK4DnqLA0jKiIiCcYnQR1oUWfvrPK4EhERkbbxRVCr61tERBKVL4L624vJ1PUtIiKJxRdBXZ2aBqhFLSIiiccXQe0sicrUdN2eJSIiCccXQQ16JrWIiCQm3wR14JnUOkctIiKJxTdBvUMtahERSUC+Ceqq1AyydiqoRUQksfgmqCvT0tX1LSIiCcc/Qa2ubxERSUC+CerAxWQKahERSSwprW1gZg8Ck4ANzrmDg8vmAAcGN+kGbHXOFYT4bAmwHagDap1zhVGqu83UohYRkUTUalADDwN3A//XsMA5d3bDvJndCmzby+fHO+c2trfAaKlKTdcQoiIiknBaDWrn3Jtmlh9qnZkZcBZwXHTLir5KdX2LiEgCivQc9VjgG+fcmhbWO+AVM3vfzKbsbUdmNsXMisysqKysLMKy9lSZmkFqfR2pdbuivm8REZGOEmlQnws8tpf1o51zhwETgCvN7JiWNnTO3eucK3TOFebm5kZY1p6+fdSlur9FRCRxtDuozSwF+AEwp6VtnHPrg68bgH8Co9r7fZGqTE0H0KAnIiKSUCJpUZ8ArHLOlYZaaWbZZpbTMA+cBCyP4PsiUpnW8ExqBbWIiCSOVoPazB4DFgIHmlmpmf04uOocmnV7m1k/M3sx+LYP8LaZfQAsBuY5516OXult823Xt4JaREQSRzhXfZ/bwvKLQyxbD0wMzq8DRkZYX9RUpqpFLSIiicdHI5MFz1HrYjIREUkgvgnqSnV9i4hIAvJdUKvrW0REEolvgrqqMajV9S0iIonDN0HdcB91pu6jFhGRBOKboP72YjIFtYiIJA7fBHV9UjLVKWm6mExERBKKb4IaGp5JrXPUIiKSOHwW1Onq+hYRkYTiq6CuSs0ga2eV12WIiIiEzVdBra5vERFJNL4K6qrUdF1MJiIiCcVXQR1oUSuoRUQkcfgwqNX1LSIiicNXQV2VmqGubxERSSi+CurKNN2eJSIiicVXQa0WtYiIJBpfBXVlagbpdbVQW+t1KSIiImHxWVAHHszBjh3eFiIiIhImXwV1wzOpFdQiIpIofBXUlQpqERFJMApqERGROOaroK7SOWoREUkwvgrqyjS1qEVEJLH4Kqh1MZmIiCQaXwW1zlGLiEii8VlQ6xy1iIgkFl8Ftbq+RUQk0fgqqNX1LSIiicZXQV2bnMLOpBQFtYiIJAxfBTUE76VWUIuISILwXVBXpmYoqEVEJGH4L6jTFNQiIpI4/BfUqRlQUeF1GSIiImHxXVBXpGfB9u1elyEiIhIW3wX19vRs2LrV6zJERETC4rugLk/Phm3bvC5DREQkLL4L6u3pWWpRi4hIwvBdUJenZwfOUdfXe12KiIhIq/wZ1M5BebnXpYiIiLTKd0G9PT07MKPz1CIikgB8F9TlGQpqERFJHK0GtZk9aGYbzGx5k2XXm9mXZrYsOE1s4bOnmNnHZrbWzKZHs/D2amxR64IyERFJAOG0qB8GTgmx/DbnXEFwerH5SjNLBv4KTACGAeea2bBIio2G7elZgRm1qEVEJAG0GtTOuTeBze3Y9yhgrXNunXNuJ/A4MLkd+4mqcrWoRUQkgURyjvoqMysOdo13D7G+P/BFk/elwWUhmdkUMysys6KysrIIyto7XUwmIiKJpL1BfQ+wP1AAfAXcGmIbC7HMtbRD59y9zrlC51xhbm5uO8tqnYJaREQSSbuC2jn3jXOuzjlXD9xHoJu7uVJgQJP3ecD69nxfNO1MSYWMDHV9i4hIQmhXUJtZ3yZvvw8sD7HZe8ABZjbIzNKAc4C57fm+qOvaVS1qERFJCCmtbWBmjwHjgF5mVgpcB4wzswICXdklwOXBbfsB9zvnJjrnas3sKuBfQDLwoHNuRYf8FG3VrZta1CIikhBaDWrn3LkhFj/QwrbrgYlN3r8I7HHrlufUohYRkQThu5HJgECLWkEtIiIJwJ9B3bWrur5FRCQh+DOo1aIWEZEE4cug/tsHm6jeuJn86fPInz7P63JERERa5Mug3p6eRUbtTlLrdnldioiIyF75MqgbxvvOqan0uBIREZG982VQNwwjuk91hceViIiI7J0vg7o8Qy1qERFJDL4M6sYWdc0OjysRERHZO58GdRYAOQpqERGJc74M6vL0LoCCWkRE4p8vg7qhRb1PtYJaRETim2+Duh5jH11MJiIicc6XQe0siYq0TPap0e1ZIiIS33wZ1BC4RUu3Z4mISLzzbVBvT8/WxWQiIhL3fBzUWbqPWkRE4p5vg7o8XV3fIiIS/3wb1NvTszXWt4iIxD3fBrUuJhMRkUTg26BuvJjMOa9LERERaZGPgzqLFFdP1q5qr0sRERFpkW+DWuN9i4hIIvBtUGu8bxERSQS+Dery4DOpdUGZiIjEM98G9fZgUGu8bxERiWe+DeryjIagVotaRETil3+DurHrW+eoRUQkfvk2qBsvJlNQi4hIHPNtUFenpLMzKUUtahERiWu+DWrMAk/Q0u1ZIiISx/wb1Gi8bxERiX++DurG8b5FRETilM+DOksXk4mISFzzdVCXp3dRi1pEROKar4NaF5OJiEi883VQl6dnk7NTF5OJiEj88nVQb0/PpsvOKqit9boUERGRkHwf1ACUl3tbiIiISAt8HdQND+Zg2zZvCxEREWmBr4O6Ybxvtm71thAREZEW+Dqoy9O7BGbUohYRkTjValCb2YNmtsHMljdZdouZrTKzYjP7p5l1a+GzJWb2oZktM7OiaBYeDeUNLWoFtYiIxKlwWtQPA6c0W/YqcLBzbgSwGvjNXj4/3jlX4JwrbF+JHac8I9iiVte3iIjEqVaD2jn3JrC52bJXnHMN9zQtAvI6oLYOt10tahERiXPROEd9KfBSC+sc8IqZvW9mU6LwXVFVkaaLyUREJL6lRPJhM/sdUAvMbmGT0c659WbWG3jVzFYFW+ih9jUFmAKw3377RVJW2GqTU9iRmkG2WtQiIhKn2t2iNrOLgEnAec45F2ob59z64OsG4J/AqJb255y71zlX6JwrzM3NbW9ZbbY9PUtd3yIiErfaFdRmdgrwa+B051zIwbLNLNvMchrmgZOA5aG29VJ5ehd1fYuISNwK5/asx4CFwIFmVmpmPwbuBnIIdGcvM7NZwW37mdmLwY/2Ad42sw+AxcA859zLHfJTREAtahERiWetnqN2zp0bYvEDLWy7HpgYnF8HjIyouhgoz8hWi1pEROKWr0cmA9iakQObNnldhoiISEi+D+rNWV2hrMzrMkRERELyfVBvyuoKFRVQXe11KSIiInvwfVBvztwnMKNWtYiIxCEFdVbXwIyCWkRE4pDvg3qTglpEROKY74NaLWoREYlnvg9qtahFRCSe+T6oy9OzISVFQS0iInHJ90GNGfTqpaAWEZG4pKAGyM1VUIuISFxSUIOCWkRE4paCGhTUIiIStxTUoKAWEZG4paCGQFBv3Qq7dnldiYiIyG4U1BAIaoCNG72tQ0REpBkFNXwb1Or+FhGROKOgBgW1iIjELQU1KKhFRCRuKahBQS0iInFLQQ3Qo0dgKFEFtYiIxBkFNUByMvTsqaAWEZG4o6BuoAdziIhIHFJQN9DoZCIiEocU1A0U1CIiEocU1A0U1CIiEocU1A1yc2HzZqir87oSERGRRgrqBrm5UF8PW7Z4XYmIiEgjBXUDDXoiIiJxSEHdQEEtIiJxSEHdQEEtIiJxSEHdQEEtIiJxSEHdoFevwKuCWkRE4oiCukFaGnTtqqAWEZG4oqBuSoOeiIhInFFQN6WgFhGROKOgbkpBLSIicUZB3ZSCWkRE4oyCuqncXNi4EZzzuhIRERFAQb273FzYtQu2bfO6EhEREUBBvTsNeiIiInFGQd2UglpEROKMgropBbWIiMSZsILazB40sw1mtrzJsh5m9qqZrQm+dm/hsxcFt1ljZhdFq/AOoaAWEZE4E26L+mHglGbLpgOvO+cOAF4Pvt+NmfUArgOOAEYB17UU6HFBQS0iInEmrKB2zr0JbG62eDLwSHD+EeB7IT56MvCqc26zc24L8Cp7Bn78yMyE7GwFtYiIxI1IzlH3cc59BRB87R1im/7AF03elwaX7cHMpphZkZkVlXkZlBr0RERE4khHX0xmIZaFHE3EOXevc67QOVeY29AF7QUFtYiIxJFIgvobM+sLEHzdEGKbUmBAk/d5wPoIvrPjKahFRCSORBLUc4GGq7gvAp4Lsc2/gJPMrHvwIrKTgsvil4JaRETiSLi3Zz0GLAQONLNSM/sxMAM40czWACcG32NmhWZ2P4BzbjPwB+C94HRDcFn86tsXvv4a6uu9rkRERISUcDZyzp3bwqrjQ2xbBFzW5P2DwIPtqs4LeXmB8b7LyqBPH6+rERERn9PIZM3l5QVeS0u9rUNERAQF9Z4U1CIiEkcU1M0pqEVEJI4oqJvLzYXUVAW1iIjEBQV1c0lJ0L+/glpEROKCgjqUvDwFtYiIxAUFdSgKahERiRMK6lAagtqFHJZcREQkZhTUoeTlQXU1bI7vQdRERKTzU1CH0j/4JE51f4uIiMcU1KHoXmoREYkTCupQGoL6yy+9rUNERHwvrIdydHb50+ft9i7HE9gAABuTSURBVL7kxpMD91OrRS0iIh5TizqUlJTA4y4V1CIi4jEFdUt0L7WIiMQBBXVLFNQiIhIHFNQtUVCLiEgcUFC3JC8Ptm+H8nKvKxERER9TULdE91KLiEgcUFC3REEtIiJxQEHdEgW1iIjEAQV1S/r1C7wqqEVExEMK6pakpUGfPgpqERHxlIYQ3ZvgLVp7DDE641SPChIREb9Ri3pvdC+1iIh4TEG9N/37K6hFRMRTCuq9ycuDLVvI3FntdSUiIuJTCuq9Cd6itW/FJo8LERERv1JQ701DUG/f6HEhIiLiVwrqvWkMarWoRUTEG7o9K4SG27EydlWzCuirFrWIiHhELeq9qE7NYEtGjlrUIiLiGQV1K77O6akWtYiIeEZB3YqvcnrpYjIREfGMgroVX+f0Ute3iIh4RkHdivX79CK3civpu2q8LkVERHxIQd2Kz7v1BWC/rV97XImIiPiRgroVJd0DQZ2/9SuPKxERET9SULeipHs/API3r/e4EhER8SMFdSvKM7qwOXMf8rcqqEVEJPYU1GEo6d6X/C0KahERiT0FdRhKuvdj4BadoxYRkdhrd1Cb2YFmtqzJVG5m1zTbZpyZbWuyzbWRlxx7Jd370a98I+m1O70uRUREfKbdD+Vwzn0MFACYWTLwJfDPEJu+5Zyb1N7viQcl3fuShGPA1q9Z22s/r8sREREfiVbX9/HAJ865z6K0v7jScOX3IJ2nFhGRGItWUJ8DPNbCuqPM7AMze8nMhre0AzObYmZFZlZUVlYWpbKioyGoByqoRUQkxiIOajNLA04Hngyxegkw0Dk3ErgLeLal/Tjn7nXOFTrnCnNzcyMtK6rKM7qwJSNHLWoREYm5aLSoJwBLnHPfNF/hnCt3zlUE518EUs2sVxS+M+Z05beIiHghGkF9Li10e5vZvmZmwflRwe9LyEdR6V5qERHxQkRBbWZZwInAM02WTTWzqcG3ZwDLzewD4E7gHOeci+Q7vaJbtERExAvtvj0LwDlXCfRstmxWk/m7gbsj+Y540fQWLRERkVjRyGRhanw4h85Ti4hIDCmow/RtUH/pcSUiIuInCuowbcvMYUtGjlrUIiISUwrqNvhMV36LiEiMKajb4NPu/dSiFhGRmFJQt8Fn3fvSr7wMqqu9LkVERHxCQd0Gn3bvRxIOPv3U61JERMQnFNRt8Fnwym/WrPG2EBER8Q0FdRt82hDUa9d6W4iIiPiGgroNtmXmsDWji1rUIiISMwrqNirp3k8tahERiRkFdRuVdO+rFrWIiMSMgrqNSrr3g88/h5oar0sREREfUFC30boeeeAcrF7tdSkiIuIDCuo2+ig3PzDzwQee1iEiIv6goG6jdT3zIC1NQS0iIjGhoG6juqRkGD5cQS0iIjGhoG6PkSMV1CIiEhMK6vYYORI2bIBvvvG6EhER6eQU1O0xcmTgVa1qERHpYArq9lBQi4hIjCio26NHD8jLU1CLiEiHU1C3ly4oExGRGFBQt9fIkbBqlYYSFRGRDqWgbq+RI6G2Flau9LoSERHpxBTU7aULykREJAYU1O31ne9AZqaCWkREOpSCur2Sk+HggxXUIiLSoRTUkRg5EoqLA4+9FBER6QAK6kiMHAmbNsH69V5XIiIinZSCOhK6oExERDqYgjoSI0YEXhXUIiLSQRTUkejaFfLzFdQiItJhFNSR0lCiIiLSgRTUkRo5ElavhqoqrysREZFOSEEdqZEjob4+cJuWiIhIlCmoI3XkkYHXd9/1tg4REemUFNSR6tcvcEHZO+94XYmIiHRCKV4XkIjyp8/b7X3JmDHw2muBEcrMPKpKREQ6I7Woo2H0aPj6a1i3zutKRESkk1FQR8Po0YFXdX+LiEiUKaijYfjwwOAnCmoREYmyiIPazErM7EMzW2ZmRSHWm5ndaWZrzazYzA6L9DvjTlISHH00vP2215WIiEgnE60W9XjnXIFzrjDEugnAAcFpCnBPlL4zvowZAytXwubNXlciIiKdSCy6vicD/+cCFgHdzKxvDL43thrOUy9c6G0dIiLSqUQjqB3wipm9b2ZTQqzvD3zR5H1pcNluzGyKmRWZWVFZWVkUyoqxww+HlBR1f4uISFRFI6hHO+cOI9DFfaWZHdNsfagbi90eC5y71zlX6JwrzM3NjUJZMZaVBd/9ri4oExGRqIo4qJ1z64OvG4B/AqOabVIKDGjyPg9YH+n3xqXRo+G996CmxutKRESkk4goqM0s28xyGuaBk4DlzTabC1wYvPr7SGCbc+6rSL43bo0eDdXVsGSJ15WIiEgnEekQon2Af1pg2MwU4B/OuZfNbCqAc24W8CIwEVgLVAKXRPid8avpwCdHHeVtLSIi0ilEFNTOuXXAyBDLZzWZd8CVkXxPvGs69vf87n0Z9M478ItfeFiRiIh0FhqZLMqK+g8PtKjdHtfLiYiItJmCOsqK8g6CsjJYs8brUkREpBNQUEfZe3nDAzNvvOFtISIi0ikoqKNsXY/+MGgQvPSS16WIiEgnoKCONjOYOBFeey1wq5aIiEgEFNQdYeJEqKyEN9/0uhIREUlwCuqOMG4cZGTAiy96XYmIiCQ4BXVHyMqC445TUIuISMQU1B1l4sTALVq6TUtERCKgoO4oEyYEXtWqFhGRCCioO8rgwTB0qIJaREQioqDuSBMnwoIFsGOH15WIiEiCUlB3pIkTYedOjVImIiLtpqDuSGPGQJcu6v4WEZF2U1B3pPR0OOEEmDdPT9MSEZF2ieh51BKGiRPh2WdhxQo4+GBg9+dXA5TMONWLykREJAGoRd3RJk4MvD77rLd1iIhIQlJQd7T+/WHsWJg9W93fIiLSZgrqWDj/fFi1CpYu9boSERFJMArqWDjjDEhNhb//3etKREQkwSioY6FHDzj1VHjsMair87oaERFJIArqWDn/fPj6aw1+IiIibaKgjpVTT4WuXQMXlYmIiIRJQR0rGRmBc9VPP03GrmqvqxERkQShoI6B/OnzyJ8+j3N37A8VFZywdrHXJYmISIJQUMfQov0OZn1OL763Yr7XpYiISILQEKIdoPkQoQ2cJTF32LH8+L1n6V65jS1ZXWNcmYiIJBq1qGPs2WHjSK2v47SP3vS6FBERSQAK6hhb1XsQxft+hwuXzMNcvdfliIhInFNQe+DBwsl8Z3MpYz/VkKIiIrJ3CmoPzBs6hg3Z3bm0aK7XpYiISJxTUHtgV3Iq/3fYqYz79H323/iF1+WIiEgcU1B75B8FE6hJTuWS99WqFhGRlimoPbI5qyvPDhvHD5e/AZs3e12OiIjEKQW1hx4qPJ3M2hq47z6vSxERkTiloPbQqt6DeGfgCLj7bti1y+tyREQkDimoPfbQdydDaSk884zXpYiISBxSUHvsjf0L4YAD4KaboF4DoIiIyO4U1B6rT0rm6qHfgw8/5KeTf9X4pC0RERFQUMeF5w8ay8reg/jZ27NJrdO5ahER+ZaCOg44S+LPx1zIwK1fc3bxq16XIyIicURBHScWDC5kcd4wrn7nMTJ2VXtdjoiIxIl2B7WZDTCz+Wb2kZmtMLOfhthmnJltM7NlwenayMrtxMz487EX0XvHFi55/3mvqxERkTiREsFna4GfO+eWmFkO8L6ZveqcW9lsu7ecc5Mi+B7fKMobzuv7H87URU/Bltuge3evSxIREY+1u0XtnPvKObckOL8d+AjoH63C/GrmMRfQtWYH90yY0ngFuK4CFxHxr6icozazfOBQ4D8hVh9lZh+Y2UtmNjwa39eZfdR7ME8dfDw/fu9ZDij7zOtyRETEYxEHtZl1AZ4GrnHOlTdbvQQY6JwbCdwFPLuX/UwxsyIzKyorK4u0rIT2x/GXUpGexZ/+dTfmNAiKiIifRRTUZpZKIKRnO+f2GAPTOVfunKsIzr8IpJpZr1D7cs7d65wrdM4V5ubmRlJWwtuc1ZUbj/sxhV9+xHnLXva6HBER8VAkV30b8ADwkXPuLy1ss29wO8xsVPD7NrX3O/3kmeHH8dbAAn614GH6bN/odTkiIuKRSFrUo4ELgOOa3H410cymmtnU4DZnAMvN7APgTuAc55yLsGZ/MON3J19JWn0t1792r9fViIiIR9p9e5Zz7m3AWtnmbuDu9n6H333evS+3j/4R0//9MDz7LHzve16XJCIiMaaRyeLc/Yd/j5W9B8HUqfD1116XIyIiMaagjnO1ySn8dNIvoLwcfvQjqKvzuiQREYkhBXUCWJM7EO65B+bPh9//3utyREQkhiIZQlRi6aKL4N//hhtvhDFj4KSTItpd89HOSmacGtH+RESkY6hFnUjuvhuGD4fzzoMvv/S6GhERiQG1qBNJVhY8+SQUFsJZZ8Hrr0NGBqAWsohIZ6UWdaIZOhQeegjefRcuuADqNcSoiEhnpqBORGeeCbfeCk89BdOmgcaQERHptNT1nah+9jMoLYXbboO8PGCY1xWJiEgHUIs6kc2cCWefDb/6Fd9bMd/rakREpAOoRZ3IkpLgkUdgwwZmzruN2qRkXjjoGK+rEhGRKFKLOtGlp8Nzz/F+3jDueH4mZxa/6nVFIiISRQrqziAnh4vOvJ53Bo7klpfu4IIlL3hdkYiIRImCupOoTs3gsh9eyysHHMkfXp0FN9+sq8FFRDoBnaPuRHampPJfk6dz67zbmDx9Onz6Kdx5J6SlRbxvDagiIuINtag7mdrkFKZN+hn85jfwt7/BCSdAWZnXZYmISDspqDuh+qRk+OMf4R//gPfeg8MPhw8+8LosERFpBwV1Z3buufDWW1BbC0cfDQ88oPPWIiIJRkHdSeVPnxeYnvqGw7/3J97J/Q5cdhkvDx1NwU8f87o8EREJk4LaB8q69OD8s2/kxvGXMv6T9/jXg1cx9tMlXpclIiJhUFD7hLMk7h/1A7534V/Ylt6FR5+4lptfvIPuldu8Lk1ERPZCt2cliOa3R7XXR70Hc9pFtzHt7dn8uOg5TlqziD+NuwTqJwSGJBURkbiiv8w+VJOazozxl3LqxXewutd+/PnlO2HsWCgq8ro0ERFpRi1qH1udm8/ZP5rBD5e/wW8WPEivww9n3oGjuXXsBazrmed1eSIiglrUYsbThxzPsVPu4/bR53Lsp0t45YH/4o8v30W/8g1eVyci4ntqUQsAO9KzuH3MeTx66KlctXAO5y19iTM/fI25w45l1qgfsiZ3YJv2pyFHRUSiQy1q2c2m7G78/oTLGXf5vTx66KlM+PgdXn3wSu5/6vcc+XmxBkwREYkxtaglpPX79OaGE6Zw5+hzuOj9F7hoyQs8/thvWdNzALMLJsDW0dCtm9dlioh0egpq2autmftwx5gfMeuIH3Laqrc4f+mLXP/6vdD/73DmmXDBBTBuHCQne11qwtPpAhEJRUEtYalJTeepQ07gqUNOYPjXa5mXtgLmzIFHHoH+/eFHPwpMI0eCmdflioh0GjpHLW22Yt/vkN9jMgde9hBXnv5rXsvox65b/wKHHspnPfrxtyN+yGFffoS5eq9LFRFJeGpRS7vVpKYz76CxzDtoLD0qt3HS6oWcvGYhlxTN5fLFz1CW3Y1/D/ouCwZ/F7YcDd27e12yiEjCUVBLVGzO6srjBafweMEp5NTs4Li1izlh7WJOWPsfzlj+OvSaCUccAccdB+PHBx67mZnpddkiInFPQS1Rtz09m+eGj+e54eNJrq9j5PrVPDNoG7z6KsyYATfdBGlpgeAeMyYQ2kcdBT17el26iEjcUVBLh6pLSmZJ3kHwh1PhD3+A8nJ4+2144w3497/hllugthaAT3rk8UHfA/hw3+9w3Q0XQ0EBZGW1+7t1FbWIdAYKaomtffaBiRMDE0BlJRQVcfO1D3LYl6s4+rMP+MGK+fD6fYGneQ0ZEgjskSMD08EHQ16eriwXEd9QUIu3srLgmGO458jtjYt6b9/E4lO6B57mtWwZLFwIjz/+7WdycmDYsMA0dGggzA88EAYPhvR0D34IEZGOo6CWmGhLN/SGnJ5w2qlw2mmNy0ZcM4ehZSUcsPHzwPT15wxZ/k9yd2z99oNJSbDffrD//rD//ly+eiefd92XL7v2prRrn8Dwpx62xNUVvycdE5HWKaglIZRndGHxgINZPODg3Zbn1Oxg0OYvGbT5S+44NAvWroV16+CZZ/jNxo277+SBnwSCfMCAQPf5gAGBqV+/wKAt/foFLmhL0vACIhI/FNSS0LanZ1PcdwjFfYdwxw27t8YOmfYE/bdtIG/bBvK2fcP1I7Lh88/hiy9g+XLqv/qaJHZ/yMjOpBQ2ZnejLLs7IwsPhH33hT59oHdv/t9rpWzK6srmrK5sztyHxX85e69d7c1biyIi7aGglk5re3o2q3oPYlXvQQBc36xb9cBfPkvvii30qdhE74rN9KnYTJ+KTeRWbCV3xxYoLQ2cJy8rg7o67mr+Bf97UeB8ec+eu0/du0P37ly2+Cu2ZXShPKML5eldKM/Ipjw9m+3pWWxPz47NQRCRhKegFk940drc4zuTU/mya2++7Np7r58zV0/X6gp67thGbuUWulVtp2flNm4a2xc2boRNm76dPvkEtmyBrVv57/pWhlD9a1bgKviuXSEnh9llO6lIy6IiLZPKtExmLXgoMJ+ayY60DG6+6GjIzoYuXQIX4WVnfztlZUFmJvm/e3mvX9n8HHCinyMOp/5E/xlFFNQirXCWxNbMfdiauQ+fMKBx+U3/vZc/+PX1HPLzp+haXcE+1TvYpybw2rW6gi47K8mpqeRnR+wL27YF7i0vLyftq88ZuOMrsndVk7Wzii47q0iv2/XtPl/eo02/h1UpaVSlpFOdkkZVajo1TV6rU9JgzYOBEeEyMiAzk98sXk9Ncio1KWnUpKTC3Z8G1qWnB6am82lpe843XZaaGph065xIVEUU1GZ2CnAHkAzc75yb0Wx9OvB/wHeBTcDZzrmSSL5TJCEkJbE9PTvQxd019CY/a9ayOzNEL0NKXS1Zu6rJ2lnNoquPgIoK2LFjz6myEiorefilD8ncVU3Grp1k1O4ks7aGjF01pNfupEdlOXz8MVRVQXU1VFdzwfYdZOza+e25+vlR+NlTUwMB3vDaMN90eXB67ItydiWlUJuUTG1yCqx96Nv1KSnfvrYwP3XRWmqTkqhLSqbOkmDWF4H1ycmN25760QfUJiVTn5QU2Obl5MD65lPD50JNSUl7Xx7qVf9gkShpd1CbWTLwV+BEoBR4z8zmOudWNtnsx8AW59x3zOwc4Gbg7EgKFokXsehSrU1OoTw5cJ6boUNb3X5GVTtOKThHSn0d6bU7WfG78YEgr6mBmhom3fIaabW7SK/bRWrdLtLqaoOvgflbJh/UuC27dnH7i8tJraslrW4XKfV1gde6OlLra0mtqyV1Zy0pVbWk1teSUl9NSn09XWorSamvI6W+jtUL1pNSX0tqfR3J9XX0y04JjFy3axfU1lJTvZPk+jpSgk9mm978Z3ltzx/vr80XPN32Q9QuZoHADhXiTaYNlbXUm1FvSdSbkdcje49tdpua7jcpiaWl5Y2fdwajBvfabbs3P9kMgDNr3O74g/rsub8mr3OLv8JhOIN6S+IHh+XtuW3TqfmycLYJNbWwzcxXVwfrCfzj51cTDtpzu4Zj3nw+1LS39c3XBd//4qniwHEMHpe/3PITOOigmPynFEmLehSw1jm3DsDMHgcmA02DejJwfXD+KeBuMzPn3O6X2oqId8yoTU4JtGh79dpt1fJ9S/b60Vsu3/0fJ7eXR/fag+b/+Dmw4R9HzpHs6kmpq20M+WRXT3J9Hcn19aS4usb55KbrXD1zpx4JdXV7TrW1UFfH5Q8vbtw+ydVzx1kjQ29fVwf19dz4/PLg9vWYc/zyxAO+Xe9c43aNrw1T8P1rC0swV09S8Gc647D+u2/j3J7vG5bV1VG+sQ5zjqTgPkhK2u17sndWkeQchsOCr3xZ+21tTfdXXw/OcXBZBYYLfM452PHpHtvsMTVfHmq7lj7bfGrmF80XvBnV/8zCMrP5ghMHxiyorb2ZaWZnAKc45y4Lvr8AOMI5d1WTbZYHtykNvv8kuM3GEPubAkwJvj0Q+LhdhXUevYA9jpNEnY5zbOg4x46OdWxE+zgPdM7lhloRSYs61AmY5qkfzjaBhc7dC9wbQT2dipkVOecKva6js9Nxjg0d59jRsY6NWB7nSIZgKoUml8BCHrC+pW3MLIXAZTWbI/hOERERX4kkqN8DDjCzQWaWBpwDzG22zVzgouD8GcAbOj8tIiISvnZ3fTvnas3sKuBfBG7PetA5t8LMbgCKnHNzgQeAR81sLYGW9DnRKNondBogNnScY0PHOXZ0rGMjZse53ReTiYiISMfTY4JERETimIJaREQkjimoPWZmp5jZx2a21sz2GGjJzH5mZivNrNjMXjezgV7UmehaO85NtjvDzJyZ6faWdgjnOJvZWcH/pleY2T9iXWNnEMbfjf3MbL6ZLQ3+7ZjoRZ2JzsweNLMNwTFBQq03M7sz+HsoNrPDOqQQ55wmjyYCF+F9AgwG0oAPgGHNthkPZAXnrwDmeF13ok3hHOfgdjkExjxaBBR6XXeiTWH+93wAsBToHnzf2+u6E20K8zjfC1wRnB8GlHhddyJOwDHAYcDyFtZPBF4iMGbIkcB/OqIOtai91TgMq3NuJ9AwDGsj59x851xl8O0iAverS9u0epyD/gD8GaiOZXGdSDjH+SfAX51zWwCccxtiXGNnEM5xdsA+wfmu7DnGhYTBOfcmex/7YzLwfy5gEdDNzPpGuw4Ftbf6A180eV8aXNaSHxP415u0TavH2cwOBQY4516IZWGdTDj/PQ8BhpjZO2a2KPgEPmmbcI7z9cD5ZlYKvAj8v9iU5jtt/RveLnoetbfCHmLVzM4HCoFjO7Sizmmvx9nMkoDbgItjVVAnFc5/zykEur/HEegdesvMDnbObe3g2jqTcI7zucDDzrlbzewoAuNZHOxc8LFjEi1h/w2PhFrU3gpnGFbM7ATgd8DpzrmaGNXWmbR2nHOAg4EFZlZC4FzTXF1Q1mbhDiv8nHNul3PuUwIP3zkgRvV1FuEc5x8DTwA45xYCGQQeIiHRFdbf8EgpqL3V6jCswS7ZvxEIaZ3Pa5+9Hmfn3DbnXC/nXL5zLp/AtQCnO+eKvCk3YYUzrPCzBC6QxMx6EegKXxfTKhNfOMf5c+B4ADM7iEBQl8W0Sn+YC1wYvPr7SGCbc+6raH+Jur495MIbhvUWoAvwpAUeYv65c+50z4pOQGEeZ4lQmMf5X8BJZrYSqAN+6Zzb5F3ViSfM4/xz4D4zm0agK/ZiF7xMWcJnZo8ROE3TK3i+/zogFcA5N4vA+f+JwFqgErikQ+rQ705ERCR+qetbREQkjimoRURE4piCWkREJI4pqEVEROKYglpERCSOKahFRETimIJaREQkjv1/g1FVMvevrNsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#define the inversed function of pdf.\n",
    "def inversed_func(f):\n",
    "    return (100 - 99*f)**(-0.5)\n",
    "\n",
    "#check sin(theta/2)\n",
    "x = inverse_sampling(inversed_func, 1000)\n",
    "ax.cla()\n",
    "ax.plot(x0, px0, 'r-', label='Theory')\n",
    "ax.hist(x, bins=100, range=(0.1, 1), density=True, label='Inversed funtion Sampling Result')\n",
    "ax.legend()\n",
    "ax.set_title(r'$\\sin\\theta/2$')\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 模拟散射过程\n",
    "\n",
    "有了抽样散射角的函数，就可以模拟散射了。注意到每个粒子发生散射的次数服从泊松分布，而我们一开始就用泊松抽样出了每个粒子的散射次数并储存在times数组里。\n",
    "\n",
    "接下来定义一个完整的进行n次散射的过程得到最后的散射角的函数，然后依次把散射次数作为参数传进去，就可以得到最终散射角的数组了。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define a pipline for simulating a whole scattering procedure and returning the final theta\n",
    "def scatter(nscatter, func, range_, mode):\n",
    "    r'''\n",
    "    Get the scattering angle for nscatter times of scattering.\n",
    "    Args: 4\n",
    "        nscatter: an integer, how many times of scattering do you need, got this from Poisson dis..\n",
    "        func: a function, the pdf. to be sampled. or the inversed function.\n",
    "        range_: a list, [[xmin, xmax], [ymin, ymax]]. or None (for inverse mode).\n",
    "        mode: a string, 'mc' or 'inverse'.\n",
    "    Return: 1\n",
    "        theta: a float, the final scattering angle after nscatter times of scattering.\n",
    "    '''\n",
    "    assert(mode in ['mc', 'inverse'])\n",
    "    if mode == 'mc':\n",
    "        x_arr = sampling(func, range_=range_, n=nscatter)\n",
    "    elif mode == 'inverse':\n",
    "        x_arr = inverse_sampling(func, n=nscatter)\n",
    "    else:\n",
    "        raise Exception('invalid mode {}'.format(mode))\n",
    "    sign_arr = np.random.rand(nscatter)\n",
    "    sign_arr = np.where(sign_arr>0.5, 1, -1)\n",
    "    theta_arr = np.arcsin(x_arr) * 2 * sign_arr\n",
    "    theta = theta_arr.sum()\n",
    "    return theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████| 10000/10000 [00:00<00:00, 11954.69it/s]\n"
     ]
    }
   ],
   "source": [
    "#simulate\n",
    "theta_arr1 = np.zeros_like(times, dtype=np.float32)\n",
    "theta_arr2 = np.zeros_like(times, dtype=np.float32)\n",
    "for idx in tqdm(range(n)):\n",
    "    theta_arr1[idx] = scatter(times[idx], func, range_, 'mc')\n",
    "    theta_arr2[idx] = scatter(times[idx], inversed_func, None, 'inverse')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAHiCAYAAADF4pQuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dfbRddX3v+/dHELTikSDRgwFNtLEVb28jIwM59d7jAxYQTw32SBuO1WhxxAdo7TjtHSdoz9FjS4tnnJZRR60WKxWtEqgP11jS0ogwbG+LElsEAkW2QE1MCrEB1DKkAt/7x/ptOtlZe++1k/00d96vMdZYc37nb871++25ks+ec809V6oKSZK0+D1hoTsgSZJGY2hLktQThrYkST1haEuS1BOGtiRJPWFoS5LUE4a2pMck+X6S5y50PyQNZ2hLekxVHVVVd850vSQvS1JJPjuh/pOtfl2nliS/nOSWJP+SZFeSP03yE7MwBGlJM7QlzZa9wE8leXqntgH4xoR2vwe8E/hl4Bjg+cD/C7x6Pjop9ZmhLfVEkruTXJDk1iT3JfnjJE9qy77QTm2PPx5N8qYh2/iLJOdPqH09yc+26Uryo2361Un+Psl3k+xM8t5puvivDMJ3fVv/MODngE92Xms1cB5wTlV9qaoeqqoHq+qTVXXRAf5opEOGoS31y+uB04HnMThC/XWAqvqZdmr7KOB1wD8B1wxZ/1PAOeMzSU4EngNcNaTtvwBvBI5mcBT89iRnTdO/j7d1aP3cAezuLD8V2FVVX51mO5KGMLSlfvn9qtpZVfuAC+kEMECS5zMIzp+vqp1D1v8csCbJc9r864HPVtVDExtW1XVVdXNVPVpVNwGXAy+dqnNV9TfAMUl+jEF4f3xCk6cDe6YdpaShDG2pX7pB/I/As8ZnkjwN+Dzw36vqr4atXFXfY3BUvb6V1tM5fd2V5MVJrk2yN8kDwNuAY0fo4yeA84GXM/gloeufgeNG2IakIQxtqV9O6Ew/m3bqOckTGJz6vraq/nCabVwOnJPkPwBPBq6dpN2ngC3ACVX1NODDQEbo4yeAdwBbq+rBCcuuAY5PsnaE7UiawNCW+uW8JMcnOQZ4F3BFq18IPIXBVdnT2crgc+z3AVdU1aOTtHsqsK+qfpDkZOC/jNLBqrqLwWn0dw9ZdgfwB8Dl7c/EjkjypCTrk2waZfvSoczQlvrlU8BfAne2x2+2+jnAKcB9nSvIXz9sA+3z688Cr2zbm8w7gPcl+R7wP4ArR+1kVf11Ve2eZPEvA78PfBC4H/gm8FrgC6NuXzpUpaoWug+SRpDkbuAtVfXFhe6LpIXhkbYkST1haEuS1BOeHpckqSc80pYkqScMbUmSeuLwhe7AVI499thauXLlQndDkqR587Wvfe07VbV82LJFHdorV65k+/btC90NSZLmTZJ/nGyZp8clSeoJQ1uSpJ4wtCVJ6glDW5KknjC0JUnqCUNbkqSeMLQlSeqJRf132pIkLTYrN131uPm7L3r1vL22R9qSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1xLShneRJSb6a5OtJdiT5n62+KslXktyR5IokR7T6kW1+rC1f2dnWBa1+e5LT52pQkiQtRaMcaT8EvKKqfhJYA5yR5BTg/cDFVbUauA84t7U/F7ivqn4UuLi1I8mJwHrghcAZwB8kOWw2ByNJ0lI2bWjXwPfb7BPbo4BXAJ9u9cuAs9r0ujZPW35qkrT65qp6qKruAsaAk2dlFJIkHQJG+kw7yWFJbgTuBbYB3wTur6qHW5NdwIo2vQLYCdCWPwA8vVsfso4kSZrGSKFdVY9U1RrgeAZHxy8Y1qw9Z5Jlk9UfJ8nGJNuTbN+7d+8o3ZMk6ZAwo6vHq+p+4DrgFODoJOO3QT0e2N2mdwEnALTlTwP2detD1um+xiVVtbaq1i5fvnwm3ZMkaUkb5erx5UmObtNPBl4J3AZcC7yuNdsAfL5Nb2nztOVfqqpq9fXt6vJVwGrgq7M1EEmSlrpRvjDkOOCydqX3E4Arq+rPktwKbE7ym8DfAx9t7T8KfCLJGIMj7PUAVbUjyZXArcDDwHlV9cjsDkeSpKVr2tCuqpuAFw2p38mQq7+r6gfA2ZNs60Lgwpl3U5IkeUc0SZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1Jknpi2tBOckKSa5PclmRHkne2+nuTfDvJje1xZmedC5KMJbk9yemd+hmtNpZk09wMSZKkpenwEdo8DPxqVf1dkqcCX0uyrS27uKr+d7dxkhOB9cALgWcBX0zy/Lb4g8BPA7uAG5JsqapbZ2MgkiQtddOGdlXtAfa06e8luQ1YMcUq64DNVfUQcFeSMeDktmysqu4ESLK5tTW0JUkawYw+006yEngR8JVWOj/JTUkuTbKs1VYAOzur7Wq1yeqSJGkEI4d2kqOAzwC/UlXfBT4EPA9Yw+BI/HfGmw5ZvaaoT3ydjUm2J9m+d+/eUbsnSdKSN1JoJ3kig8D+ZFV9FqCq7qmqR6rqUeAj/Nsp8F3ACZ3Vjwd2T1F/nKq6pKrWVtXa5cuXz3Q8kiQtWaNcPR7go8BtVfW7nfpxnWavBW5p01uA9UmOTLIKWA18FbgBWJ1kVZIjGFystmV2hiFJ0tI3ytXjLwHeANyc5MZWexdwTpI1DE5x3w28FaCqdiS5ksEFZg8D51XVIwBJzgeuBg4DLq2qHbM4FkmSlrRRrh7/a4Z/Hr11inUuBC4cUt861XqSJGly3hFNkqSeMLQlSeoJQ1uSpJ4wtCVJ6glDW5KknjC0JUnqCUNbkqSeMLQlSeoJQ1uSpJ4wtCVJ6glDW5KknjC0JUnqCUNbkqSeMLQlSeoJQ1uSpJ4wtCVJ6glDW5KknjC0JUnqCUNbkqSeMLQlSeoJQ1uSpJ6YNrSTnJDk2iS3JdmR5J2tfkySbUnuaM/LWj1JPpBkLMlNSU7qbGtDa39Hkg1zNyxJkpaeUY60HwZ+tapeAJwCnJfkRGATcE1VrQauafMArwJWt8dG4EMwCHngPcCLgZOB94wHvSRJmt60oV1Ve6rq79r094DbgBXAOuCy1uwy4Kw2vQ74eA1cDxyd5DjgdGBbVe2rqvuAbcAZszoaSZKWsBl9pp1kJfAi4CvAM6tqDwyCHXhGa7YC2NlZbVerTVaXJEkjGDm0kxwFfAb4lar67lRNh9RqivrE19mYZHuS7Xv37h21e5IkLXkjhXaSJzII7E9W1Wdb+Z522pv2fG+r7wJO6Kx+PLB7ivrjVNUlVbW2qtYuX758JmORJGlJO3y6BkkCfBS4rap+t7NoC7ABuKg9f75TPz/JZgYXnT1QVXuSXA38Vufis9OAC2ZnGJIkzb6Vm65a6C48zrShDbwEeANwc5IbW+1dDML6yiTnAt8Czm7LtgJnAmPAg8CbAapqX5LfAG5o7d5XVftmZRSSJB0Cpg3tqvprhn8eDXDqkPYFnDfJti4FLp1JByVJ0oB3RJMkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJ6YN7SSXJrk3yS2d2nuTfDvJje1xZmfZBUnGktye5PRO/YxWG0uyafaHIknS0jbKkfbHgDOG1C+uqjXtsRUgyYnAeuCFbZ0/SHJYksOADwKvAk4EzmltJUnSiA6frkFVfTnJyhG3tw7YXFUPAXclGQNObsvGqupOgCSbW9tbZ9xjSZIOUQfzmfb5SW5qp8+XtdoKYGenza5Wm6y+nyQbk2xPsn3v3r0H0T1JkpaWAw3tDwHPA9YAe4DfafUMaVtT1PcvVl1SVWurau3y5csPsHuSJC09054eH6aq7hmfTvIR4M/a7C7ghE7T44HdbXqyuiRJGsEBHWknOa4z+1pg/MryLcD6JEcmWQWsBr4K3ACsTrIqyREMLlbbcuDdliTp0DPtkXaSy4GXAccm2QW8B3hZkjUMTnHfDbwVoKp2JLmSwQVmDwPnVdUjbTvnA1cDhwGXVtWOWR+NJElL2ChXj58zpPzRKdpfCFw4pL4V2Dqj3kmSpMd4RzRJknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqSemDe0klya5N8ktndoxSbYluaM9L2v1JPlAkrEkNyU5qbPOhtb+jiQb5mY4kiQtXaMcaX8MOGNCbRNwTVWtBq5p8wCvAla3x0bgQzAIeeA9wIuBk4H3jAe9JEkazbShXVVfBvZNKK8DLmvTlwFndeofr4HrgaOTHAecDmyrqn1VdR+wjf1/EZAkSVM40M+0n1lVewDa8zNafQWws9NuV6tNVpckSSOa7QvRMqRWU9T330CyMcn2JNv37t07q52TJKnPDjS072mnvWnP97b6LuCETrvjgd1T1PdTVZdU1dqqWrt8+fID7J4kSUvPgYb2FmD8CvANwOc79Te2q8hPAR5op8+vBk5LsqxdgHZaq0mSpBEdPl2DJJcDLwOOTbKLwVXgFwFXJjkX+BZwdmu+FTgTGAMeBN4MUFX7kvwGcENr976qmnhxmyRJmsK0oV1V50yy6NQhbQs4b5LtXApcOqPeSZKkx3hHNEmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJ6a9uYokSYeKlZuuWuguTMkjbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqScOKrST3J3k5iQ3Jtneasck2Zbkjva8rNWT5ANJxpLclOSk2RiAJEmHitk40n55Va2pqrVtfhNwTVWtBq5p8wCvAla3x0bgQ7Pw2pIkHTLm4vT4OuCyNn0ZcFan/vEauB44Oslxc/D6kiQtSQcb2gX8ZZKvJdnYas+sqj0A7fkZrb4C2NlZd1erSZKkERx+kOu/pKp2J3kGsC3JP0zRNkNqtV+jQfhvBHj2s599kN2TJGnpOKgj7ara3Z7vBT4HnAzcM37auz3f25rvAk7orH48sHvINi+pqrVVtXb58uUH0z1JkpaUAw7tJE9J8tTxaeA04BZgC7ChNdsAfL5NbwHe2K4iPwV4YPw0uiRJmt7BnB5/JvC5JOPb+VRV/UWSG4Ark5wLfAs4u7XfCpwJjAEPAm8+iNeWJOmgrdx01UJ3YUYOOLSr6k7gJ4fU/xk4dUi9gPMO9PUkSTpYfQvpibwjmiRJPWFoS5LUEwf7J1+SJC1afT8dPpFH2pIk9YShLUlSTxjakiT1hKEtSVJPeCGaJGnJWGoXnk3kkbYkST3hkbYkqbeW+pH1RB5pS5LUE4a2JEk94elxSVJvHGqnwyfySFuSpJ7wSFuStGgd6kfWExnakqRFw5CemqfHJUnqCY+0JUlzxiPn2WVoS5JmjSE9twxtSdKkJobw3Re9esrlmluGtiRpZIb0wpr30E5yBvB7wGHAH1XVRfPdB0k6FEx3lDxdey0+8xraSQ4DPgj8NLALuCHJlqq6dT77IUnzYaahOdP1ZxqyhnL/zfeR9snAWFXdCZBkM7AOMLQl9c58h6ahq/kO7RXAzs78LuDF89yHx5nuH8FCX3RxIKezpuvzTMc00z4c7Pb9j0mShktVzd+LJWcDp1fVW9r8G4CTq+qXOm02Ahvb7I8Bt89yN44FvjPL21wIS2Uc4FgWq6UylqUyDnAsi9Vsj+U5VbV82IL5PtLeBZzQmT8e2N1tUFWXAJfMVQeSbK+qtXO1/fmyVMYBjmWxWipjWSrjAMeyWM3nWOb7NqY3AKuTrEpyBLAe2DLPfZAkqZfm9Ui7qh5Ocj5wNYM/+bq0qnbMZx8kSeqref877araCmyd79ftmLNT7/NsqYwDHMtitVTGslTGAY5lsZq3sczrhWiSJOnA+dWckiT1xJIM7SRnJ9mR5NEkk17Rl+SMJLcnGUuyqVNfleQrSe5IckW7aG7eJTkmybbWj21Jlg1p8/IkN3YeP0hyVlv2sSR3dZatmf9RPNbPacfS2j3S6e+WTn1R7JPWl1H2y5okf9vehzcl+fnOsgXdL5O97zvLj2w/47H2M1/ZWXZBq9+e5PT57PcwI4zlvya5te2Da5I8p7Ns6HttoYwwljcl2dvp81s6yza09+MdSTbMb8/36+d047i4M4ZvJLm/s2yx7ZNLk9yb5JZJlifJB9pYb0pyUmfZ3OyTqlpyD+AFDP7G+zpg7SRtDgO+CTwXOAL4OnBiW3YlsL5Nfxh4+wKN438Bm9r0JuD907Q/BtgH/Eib/xjwuoXeHzMZC/D9SeqLYp+MOhbg+cDqNv0sYA9w9ELvl6ne95027wA+3KbXA1e06RNb+yOBVW07hy3gfhhlLC/v/Ht4+/hYpnqvLeKxvAn4/SHrHgPc2Z6Xtelli3UcE9r/EoMLkhfdPmn9+Y/AScAtkyw/E/hzIMApwFfmep8sySPtqrqtqqa7Kctjt1Stqn8FNgPrkgR4BfDp1u4y4Ky56+2U1rXXH7UfrwP+vKoenNNeHZiZjuUxi2yfwAhjqapvVNUdbXo3cC8w9GYJ82zo+35Cm+74Pg2c2vbBOmBzVT1UVXcBY217C2XasVTVtZ1/D9czuDfEYjTKfpnM6cC2qtpXVfcB24Az5qif05npOM4BLp+Xnh2AqvoygwOhyawDPl4D1wNHJzmOOdwnSzK0RzTslqorgKcD91fVwxPqC+GZVbUHoD0/Y5r269n/H8CF7bTNxUmOnItOjmjUsTwpyfYk14+f5mdx7ROY4X5JcjKDo45vdsoLtV8me98PbdN+5g8w2AejrDufZtqfcxkcFY0b9l5bKKOO5T+3982nk4zfqGox7ZeR+9I+qlgFfKlTXkz7ZBSTjXfO9klvv087yReBfz9k0bur6vOjbGJIraaoz4mpxjHD7RwH/ASDv4EfdwHwTwwC4xLgvwHvO7CejtSH2RjLs6tqd5LnAl9KcjPw3SHt5vTPHmZ5v3wC2FBVj7byvO6XiV0aUpv4s1wU/zZGMHJ/kvwCsBZ4aae833utqr45bP15MMpYvgBcXlUPJXkbg7Mhrxhx3fkyk76sBz5dVY90aotpn4xi3v+t9Da0q+qVB7mJyW6p+h0GpzgOb0cZ+91qdTZNNY4k9yQ5rqr2tP/8751iUz8HfK6qftjZ9p42+VCSPwZ+bVY6PYnZGEs7lUxV3ZnkOuBFwGeYx33SXv+gx5Lk3wFXAb/eTp2Nb3te98sE095KuNNmV5LDgacxOEU4yrrzaaT+JHklg1+2XlpVD43XJ3mvLVRAjHKL53/uzH4EeH9n3ZdNWPe6We/haGbyHlkPnNctLLJ9MorJxjtn++RQPj0+9JaqNbiK4FoGnw8DbABGOXKfC1va64/Sj/0+G2qBMv6Z8FnA0Csg58m0Y0mybPxUcZJjgZcAty6yfQKjjeUI4HMMPu/60wnLFnK/jHIr4e74Xgd8qe2DLcD6DK4uXwWsBr46T/0eZtqxJHkR8IfAa6rq3k596Htt3nq+v1HGclxn9jXAbW36auC0NqZlwGk8/ozbfBrpVtVJfozBBVp/26kttn0yii3AG9tV5KcAD7Rfyudun8z11XcL8QBey+A3nYeAe4CrW/1ZwNZOuzOBbzD4Te7dnfpzGfxnNAb8KXDkAo3j6cA1wB3t+ZhWXwv8UafdSuDbwBMmrP8l4GYGofAnwFELuE+mHQvwU62/X2/P5y62fTKDsfwC8EPgxs5jzWLYL8Pe9wxOz7+mTT+p/YzH2s/8uZ11393Wux141ULtgxmM5Yvt/4DxfbBluvfaIh7LbwM7Wp+vBX68s+4vtv01Brx5MY+jzb8XuGjCeotxn1zO4C8/fsggU84F3ga8rS0P8ME21pvp/LXSXO0T74gmSVJPHMqnxyVJ6hVDW5KknjC0JUnqCUNbkqSeMLQlSeoJQ1uSpJ4wtCVJ6glDW1qikny/3cN5puv930mm+5Y8SQvAm6tIktQTHmlLWnDt3s3+fyRNw38k0gJIcneSC5LcmuS+JH+c5Elt2Rfaqe3xx6NJ3jRkG3+R5PwJta8n+dk2XUl+tE2/OsnfJ/lukp1J3jtF316WZNeEvv5a+x7nB5Jc0enrbUn+U6ft4Um+k+SkNn9Kkr9Jcn/r28s6ba9LcmGS/w94EHhukjcluTPJ95LcleT1nfa/2F7vviRXZ/B9zNIhxdCWFs7rgdOB5wHPB34doKp+pqqOqqqjGHzL1j8x+GKSiT7F4NvdAEhyIvAcBl8HOtG/AG8EjgZeDbw9yVkz6OvPAWcAq4D/E3hTq1/e7UMbz3eq6u+SrGh9+U3gGAZfQfqZJMs77d8AbASeCuwFPsDgi0ieyuALJG5sYzsLeBfws8By4K+Y8K120qHA0JYWzu9X1c6q2gdcyOPDjyTPBz4O/HxV7Ryy/ueANZ0jztcDn63Od0aPq6rrqurmqnq0qm5iEHgvnUFfP1BVu1tfvwCsafVPAa9J8iNt/r+0Ggy+6WxrVW1tr7sN2M7gW6DGfayqdtTge9IfBh4F/o8kT66qPVW1o7V7K/DbVXVba/tbE8YuHRIMbWnhdIP4Hxl8dSwASZ7G4Hu6/3tV/dWwlavqewyOZNe30nrgk8PaJnlxkmuT7E3yAIOvFzx2Bn39p870g8BRrQ9jDL7X+WdacL+Gfwvt5wBnt1Pj9ye5H/i/gO73Qj/2M6iqfwF+vvVtT5Krkvx4Z1u/19nOPgZfi7hiBmOQes/QlhbOCZ3pZwO7AdoFWZ8Crq2qP5xmG5cD5yT5D8CTGXzP8jCfArYAJ1TV04APMwi92TB+inwdcGsLchgE8ieq6ujO4ylVdVFn3cf9+UpVXV1VP80g2P8B+EhnW2+dsK0nV9XfzNIYpF4wtKWFc16S45Mcw+Dz2ita/ULgKcA7R9jGVgZHoe8DrqiqRydp91RgX1X9IMnJDE5jz5bNwGnA2/m3o2yAP2FwBH56ksOSPKld5Hb8sI0keWaS1yR5CvAQ8H3gkbb4w8AFSV7Y2j4tydmzOAapFwxtaeF8CvhL4M72+M1WPwc4BbivcwX564dtoH1+/VnglTw+MCd6B/C+JN8D/gdw5ewMAapqD/C3DC4cu6JT38ng6PtdDC4y2wn8P0z+/84TgF9lcMZhH4PP3N/RtvU54P3A5iTfBW4BXjVbY5D6wpurSAsgyd3AW6rqiwvdF0n94ZG2JEk9YWhLktQTnh6XJKknPNKWJKknDG1Jknri8IXuwFSOPfbYWrly5UJ3Q5KkefO1r33tO1W1fNiyRR3aK1euZPv27QvdDUmS5k2Sf5xsmafHJUnqCUNbkqSeMLQlSeoJQ1uSpJ4wtCVJ6glDW5KknjC0JUnqCUNbkqSeWNQ3V5EkabFZuemqx83ffdGr5+21PdKWJKknDG1JknrC0JYkqScMbUmSesLQliSpJ6YN7SRPSvLVJF9PsiPJ/2z1VUm+kuSOJFckOaLVj2zzY235ys62Lmj125OcPleDkiRpKRrlSPsh4BVV9ZPAGuCMJKcA7wcurqrVwH3Aua39ucB9VfWjwMWtHUlOBNYDLwTOAP4gyWGzORhJkpayaUO7Br7fZp/YHgW8Avh0q18GnNWm17V52vJTk6TVN1fVQ1V1FzAGnDwro5Ak6RAw0mfaSQ5LciNwL7AN+CZwf1U93JrsAla06RXAToC2/AHg6d36kHUkSdI0RgrtqnqkqtYAxzM4On7BsGbtOZMsm6z+OEk2JtmeZPvevXtH6Z4kSYeEGV09XlX3A9cBpwBHJxm/DerxwO42vQs4AaAtfxqwr1sfsk73NS6pqrVVtXb58uUz6Z4kSUvaKFePL09ydJt+MvBK4DbgWuB1rdkG4PNtekubpy3/UlVVq69vV5evAlYDX52tgUiStNSN8oUhxwGXtSu9nwBcWVV/luRWYHOS3wT+Hvhoa/9R4BNJxhgcYa8HqKodSa4EbgUeBs6rqkdmdziSJC1d04Z2Vd0EvGhI/U6GXP1dVT8Azp5kWxcCF868m5IkyTuiSZLUE4a2JEk9YWhLktQThrYkST1haEuS1BOGtiRJPWFoS5LUE4a2JEk9YWhLktQThrYkST1haEuS1BOGtiRJPWFoS5LUE4a2JEk9YWhLktQThrYkST1haEuS1BOGtiRJPWFoS5LUE4a2JEk9YWhLktQT04Z2khOSXJvktiQ7kryz1d+b5NtJbmyPMzvrXJBkLMntSU7v1M9otbEkm+ZmSJIkLU2Hj9DmYeBXq+rvkjwV+FqSbW3ZxVX1v7uNk5wIrAdeCDwL+GKS57fFHwR+GtgF3JBkS1XdOhsDkSRpqZs2tKtqD7CnTX8vyW3AiilWWQdsrqqHgLuSjAEnt2VjVXUnQJLNra2hLUnSCGb0mXaSlcCLgK+00vlJbkpyaZJlrbYC2NlZbVerTVaXJEkjGDm0kxwFfAb4lar6LvAh4HnAGgZH4r8z3nTI6jVFfeLrbEyyPcn2vXv3jto9SZKWvJFCO8kTGQT2J6vqswBVdU9VPVJVjwIf4d9Oge8CTuisfjywe4r641TVJVW1tqrWLl++fKbjkSRpyRrl6vEAHwVuq6rf7dSP6zR7LXBLm94CrE9yZJJVwGrgq8ANwOokq5IcweBitS2zMwxJkpa+Ua4efwnwBuDmJDe22ruAc5KsYXCK+27grQBVtSPJlQwuMHsYOK+qHgFIcj5wNXAYcGlV7ZjFsUiStKSNcvX4XzP88+itU6xzIXDhkPrWqdaTJEmT845okiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2gaapagAAA9/SURBVJIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPXEtKGd5IQk1ya5LcmOJO9s9WOSbEtyR3te1upJ8oEkY0luSnJSZ1sbWvs7kmyYu2FJkrT0jHKk/TDwq1X1AuAU4LwkJwKbgGuqajVwTZsHeBWwuj02Ah+CQcgD7wFeDJwMvGc86CVJ0vSmDe2q2lNVf9emvwfcBqwA1gGXtWaXAWe16XXAx2vgeuDoJMcBpwPbqmpfVd0HbAPOmNXRSJK0hM3oM+0kK4EXAV8BnllVe2AQ7MAzWrMVwM7OartabbL6xNfYmGR7ku179+6dSfckSVrSRg7tJEcBnwF+paq+O1XTIbWaov74QtUlVbW2qtYuX7581O5JkrTkjRTaSZ7IILA/WVWfbeV72mlv2vO9rb4LOKGz+vHA7inqkiRpBIdP1yBJgI8Ct1XV73YWbQE2ABe158936ucn2czgorMHqmpPkquB3+pcfHYacMHsDEOSpNm3ctNVC92Fx5k2tIGXAG8Abk5yY6u9i0FYX5nkXOBbwNlt2VbgTGAMeBB4M0BV7UvyG8ANrd37qmrfrIxCkqRDwLShXVV/zfDPowFOHdK+gPMm2dalwKUz6aAkSRrwjmiSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk9MG9pJLk1yb5JbOrX3Jvl2khvb48zOsguSjCW5PcnpnfoZrTaWZNPsD0WSpKVtlCPtjwFnDKlfXFVr2mMrQJITgfXAC9s6f5DksCSHAR8EXgWcCJzT2kqSpBEdPl2DqvpykpUjbm8dsLmqHgLuSjIGnNyWjVXVnQBJNre2t864x5IkHaIO5jPt85Pc1E6fL2u1FcDOTptdrTZZXZIkjehAQ/tDwPOANcAe4HdaPUPa1hT1/STZmGR7ku179+49wO5JkrT0HFBoV9U9VfVIVT0KfIR/OwW+Czih0/R4YPcU9WHbvqSq1lbV2uXLlx9I9yRJWpIOKLSTHNeZfS0wfmX5FmB9kiOTrAJWA18FbgBWJ1mV5AgGF6ttOfBuS5J06Jn2QrQklwMvA45Nsgt4D/CyJGsYnOK+G3grQFXtSHIlgwvMHgbOq6pH2nbOB64GDgMuraodsz4aSZKWsFGuHj9nSPmjU7S/ELhwSH0rsHVGvZMkSY/xjmiSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9cS0oZ3k0iT3JrmlUzsmybYkd7TnZa2eJB9IMpbkpiQnddbZ0NrfkWTD3AxHkqSla5Qj7Y8BZ0yobQKuqarVwDVtHuBVwOr22Ah8CAYhD7wHeDFwMvCe8aCXJEmjmTa0q+rLwL4J5XXAZW36MuCsTv3jNXA9cHSS44DTgW1Vta+q7gO2sf8vApIkaQoH+pn2M6tqD0B7fkarrwB2dtrtarXJ6pIkaUSzfSFahtRqivr+G0g2JtmeZPvevXtntXOSJPXZgYb2Pe20N+353lbfBZzQaXc8sHuK+n6q6pKqWltVa5cvX36A3ZMkaek50NDeAoxfAb4B+Hyn/sZ2FfkpwAPt9PnVwGlJlrUL0E5rNUmSNKLDp2uQ5HLgZcCxSXYxuAr8IuDKJOcC3wLObs23AmcCY8CDwJsBqmpfkt8Abmjt3ldVEy9ukyRJU5g2tKvqnEkWnTqkbQHnTbKdS4FLZ9Q7SZL0GO+IJklSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPTHtHNEmSDhUrN1210F2YkkfakiT1hKEtSVJPGNqSJPWEoS1JUk8Y2pIk9YShLUlSTxjakiT1hKEtSVJPGNqSJPWEoS1JUk8cVGgnuTvJzUluTLK91Y5Jsi3JHe15WasnyQeSjCW5KclJszEASZIOFbNxpP3yqlpTVWvb/CbgmqpaDVzT5gFeBaxuj43Ah2bhtSVJOmTMxenxdcBlbfoy4KxO/eM1cD1wdJLj5uD1JUlakg42tAv4yyRfS7Kx1Z5ZVXsA2vMzWn0FsLOz7q5WkyRJIzjYr+Z8SVXtTvIMYFuSf5iibYbUar9Gg/DfCPDsZz/7ILsnSdLScVBH2lW1uz3fC3wOOBm4Z/y0d3u+tzXfBZzQWf14YPeQbV5SVWurau3y5csPpnuSJC0pBxzaSZ6S5Knj08BpwC3AFmBDa7YB+Hyb3gK8sV1FfgrwwPhpdEmSNL2DOT3+TOBzSca386mq+oskNwBXJjkX+BZwdmu/FTgTGAMeBN58EK8tSdIh54BDu6ruBH5ySP2fgVOH1As470BfT5KkQ93BXogmSVJvrdx01UJ3YUYMbUnSIaNvIT2R9x6XJKknPNKWJC1ZfT+ynsgjbUmSesLQliSpJwxtSZJ6wtCWJKknvBBNkrRkLLULzybySFuSpJ4wtCVJ6glPj0uSemupnw6fyCNtSZJ6wtCWJKknPD0uSVq0DrXT39PxSFuSpJ7wSFuStGh4ZD01Q1uStGAM6ZkxtCVJs8YQnluGtiRpZBND+e6LXr1APTk0GdqSpElNd+TskfX8mvfQTnIG8HvAYcAfVdVF890HSeqjmR7lHmx7LT7zGtpJDgM+CPw0sAu4IcmWqrp1PvshSYvBdKE620e5hnL/zfeR9snAWFXdCZBkM7AOMLQlLTrThdxMQ/ZgX0+a79BeAezszO8CXjyfHZjpb7aH4kUWB/vb/1z/zA7kP7bZ7tPB/oy0NLifNd9SVfP3YsnZwOlV9ZY2/wbg5Kr6pU6bjcDGNvtjwO2z3I1jge/M8jYXwlIZBziWxWqpjGWpjAMcy2I122N5TlUtH7Zgvo+0dwEndOaPB3Z3G1TVJcAlc9WBJNurau1cbX++LJVxgGNZrJbKWJbKOMCxLFbzOZb5vvf4DcDqJKuSHAGsB7bMcx8kSeqleT3SrqqHk5wPXM3gT74uraod89kHSZL6at7/TruqtgJb5/t1O+bs1Ps8WyrjAMeyWC2VsSyVcYBjWazmbSzzeiGaJEk6cH6ftiRJPbEkQzvJ2Ul2JHk0yaRX9CU5I8ntScaSbOrUVyX5SpI7klzRLpqbd0mOSbKt9WNbkmVD2rw8yY2dxw+SnNWWfSzJXZ1la+Z/FI/1c9qxtHaPdPq7pVNfFPuk9WWU/bImyd+29+FNSX6+s2xB98tk7/vO8iPbz3is/cxXdpZd0Oq3Jzl9Pvs9zAhj+a9Jbm374Jokz+ksG/peWygjjOVNSfZ2+vyWzrIN7f14R5IN89vz/fo53Tgu7ozhG0nu7yxbbPvk0iT3JrllkuVJ8oE21puSnNRZNjf7pKqW3AN4AYO/8b4OWDtJm8OAbwLPBY4Avg6c2JZdCaxv0x8G3r5A4/hfwKY2vQl4/zTtjwH2AT/S5j8GvG6h98dMxgJ8f5L6otgno44FeD6wuk0/C9gDHL3Q+2Wq932nzTuAD7fp9cAVbfrE1v5IYFXbzmELuB9GGcvLO/8e3j4+lqnea4t4LG8Cfn/IuscAd7bnZW162WIdx4T2v8TgguRFt09af/4jcBJwyyTLzwT+HAhwCvCVud4nS/JIu6puq6rpbsry2C1Vq+pfgc3AuiQBXgF8urW7DDhr7no7pXXt9Uftx+uAP6+qB+e0VwdmpmN5zCLbJzDCWKrqG1V1R5veDdwLDL1Zwjwb+r6f0KY7vk8Dp7Z9sA7YXFUPVdVdwFjb3kKZdixVdW3n38P1DO4NsRiNsl8mczqwrar2VdV9wDbgjDnq53RmOo5zgMvnpWcHoKq+zOBAaDLrgI/XwPXA0UmOYw73yZIM7RENu6XqCuDpwP1V9fCE+kJ4ZlXtAWjPz5im/Xr2/wdwYTttc3GSI+eikyMadSxPSrI9yfXjp/lZXPsEZrhfkpzM4Kjjm53yQu2Xyd73Q9u0n/kDDPbBKOvOp5n251wGR0Xjhr3XFsqoY/nP7X3z6STjN6paTPtl5L60jypWAV/qlBfTPhnFZOOds33S2+/TTvJF4N8PWfTuqvr8KJsYUqsp6nNiqnHMcDvHAT/B4G/gx10A/BODwLgE+G/A+w6spyP1YTbG8uyq2p3kucCXktwMfHdIuzn9s4dZ3i+fADZU1aOtPK/7ZWKXhtQm/iwXxb+NEYzcnyS/AKwFXtop7/deq6pvDlt/Howyli8Al1fVQ0nexuBsyCtGXHe+zKQv64FPV9Ujndpi2iejmPd/K70N7ap65UFuYrJbqn6HwSmOw9tRxn63Wp1NU40jyT1JjquqPe0//3un2NTPAZ+rqh92tr2nTT6U5I+BX5uVTk9iNsbSTiVTVXcmuQ54EfAZ5nGftNc/6LEk+XfAVcCvt1Nn49ue1/0ywbS3Eu602ZXkcOBpDE4RjrLufBqpP0leyeCXrZdW1UPj9UneawsVEKPc4vmfO7MfAd7fWfdlE9a9btZ7OJqZvEfWA+d1C4tsn4xisvHO2T45lE+PD72lag2uIriWwefDABuAUY7c58KW9vqj9GO/z4ZaoIx/JnwWMPQKyHky7ViSLBs/VZzkWOAlwK2LbJ/AaGM5Avgcg8+7/nTCsoXcL6PcSrg7vtcBX2r7YAuwPoOry1cBq4GvzlO/h5l2LEleBPwh8JqqurdTH/pem7ee72+UsRzXmX0NcFubvho4rY1pGXAajz/jNp9GulV1kh9jcIHW33Zqi22fjGIL8MZ2FfkpwAPtl/K52ydzffXdQjyA1zL4Tech4B7g6lZ/FrC10+5M4BsMfpN7d6f+XAb/GY0BfwocuUDjeDpwDXBHez6m1dcCf9RptxL4NvCECet/CbiZQSj8CXDUAu6TaccC/FTr79fb87mLbZ/MYCy/APwQuLHzWLMY9suw9z2D0/OvadNPaj/jsfYzf25n3Xe39W4HXrVQ+2AGY/li+z9gfB9sme69tojH8tvAjtbna4Ef76z7i21/jQFvXszjaPPvBS6asN5i3CeXM/jLjx8yyJRzgbcBb2vLA3ywjfVmOn+tNFf7xDuiSZLUE4fy6XFJknrF0JYkqScMbUmSesLQliSpJwxtSZJ6wtCWJKknDG1JknrC0JYkqSf+f9qvm2kbDHEyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#check\n",
    "pz1 = np.cos(theta_arr1)\n",
    "pz2 = np.cos(theta_arr2)\n",
    "fig.clf()\n",
    "ax1 = fig.add_subplot(211)\n",
    "ax1.hist(pz1, bins=100, range=(-1, 1))\n",
    "ax1.set_title('pz via MC')\n",
    "ax2 = fig.add_subplot(212)\n",
    "ax2.hist(pz2, bins=100, range=(-1, 1))\n",
    "ax2.set_title('pz via inverse')\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
